Starting replication at:  2025-06-08 09:13:32.348567 
R version 4.4.0 (2024-04-24)
Platform: aarch64-apple-darwin20
Running under: macOS Sonoma 14.4.1

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Chicago
tzcode source: internal

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] tinytable_0.3.0      haven_2.5.4          DeclareDesign_1.0.10 fabricatr_1.0.2      randomizr_1.0.0     
 [6] blockTools_0.6.4     nbpMatching_1.5.5    emmeans_1.10.3       ggh4x_0.3.0          ggpubr_0.6.0        
[11] scales_1.4.0         estimatr_1.0.4       kableExtra_1.4.0     lubridate_1.9.3      forcats_1.0.0       
[16] stringr_1.5.1        dplyr_1.1.4          purrr_1.0.4          readr_2.1.5          tidyr_1.3.1         
[21] tibble_3.2.1         ggplot2_3.5.2        tidyverse_2.0.0     

loaded via a namespace (and not attached):
 [1] sjlabelled_1.2.0   tidyselect_1.2.1   viridisLite_0.4.2  farver_2.1.2       fastmap_1.2.0      TH.data_1.1-2     
 [7] digest_0.6.35      rpart_4.1.23       estimability_1.5.1 timechange_0.3.0   lifecycle_1.0.4    cluster_2.1.6     
[13] survival_3.5-8     magrittr_2.0.3     compiler_4.4.0     rlang_1.1.6        Hmisc_5.1-3        tools_4.4.0       
[19] data.table_1.17.0  knitr_1.46         ggsignif_0.6.4     htmlwidgets_1.6.4  xml2_1.3.6         RColorBrewer_1.1-3
[25] multcomp_1.4-26    abind_1.4-5        withr_3.0.2        foreign_0.8-86     nnet_7.3-19        colorspace_2.1-0  
[31] xtable_1.8-4       MASS_7.3-60.2      insight_1.2.0      cli_3.6.5          mvtnorm_1.2-5      rmarkdown_2.27    
[37] generics_0.1.3     rstudioapi_0.16.0  tzdb_0.4.0         splines_4.4.0      base64enc_0.1-3    vctrs_0.6.5       
[43] Matrix_1.7-0       sandwich_3.1-0     jsonlite_1.8.8     carData_3.0-5      car_3.1-2          hms_1.1.3         
[49] rstatix_0.7.2      Formula_1.2-5      htmlTable_2.4.2    systemfonts_1.1.0  qualtRics_3.2.0    glue_1.8.0        
[55] codetools_0.2-20   stringi_1.8.7      gtable_0.3.6       pillar_1.10.2      htmltools_0.5.8.1  R6_2.6.1          
[61] evaluate_0.23      lattice_0.22-6     backports_1.4.1    broom_1.0.6        Rcpp_1.0.12        checkmate_2.3.1   
[67] svglite_2.1.3      coda_0.19-4.1      gridExtra_2.3      xfun_0.44          zoo_1.8-12         pkgconfig_2.0.3   

> #################################################################################
> # Replication file for:                                                         #
> # "Balancing Precision and Retention in Experimental Design"                    #
> #                                                                               #
> # Gustavo Diaz                                                                  #
> # Northwestern University                                                       #
> # gustavo.diaz@northwestern.edu                                                 #
> #                                                                               #
> # Erin L. Rossiter                                                              #
> # University of Notre Dame                                                      #
> # erossite@nd.edu                                                               #
> #                                                                               #
> # This file analyzes the handcoding of published experiments to produce the     #
> # statistic reported in Section 1 and the statistics and table in Appendix D.   #
> #################################################################################
> 
> # Data -----
> df <- read.csv("./data/raw_data/handcoding.csv")

> # Appendix D.1 text -----
> # number of unique articles
> df %>%
+   distinct(doi) %>%
+   nrow()
[1] 221

> # number of unique articles handcoded
> df %>%
+   filter(in_sample == 1) %>%
+   distinct(doi) %>%
+   nrow()
[1] 122

> # number of unique experiments handcoded
> df %>%
+   filter(in_sample == 1) %>%
+   distinct(doi, experiment_id) %>%
+   nrow()
[1] 217

> # percent that collected pre-treatment covariates
> # among non-prepost and non-blocking experiments
> df %>%
+   filter(in_sample == 1 & block == 0 & prepost == 0) %>%
+   summarize(collect_pt = mean(pretreatment_covars))
  collect_pt
1  0.9243243

> # Appendix Table D1 -----
> design_stats <- df %>%
+   filter(in_sample == 1) %>%
+   group_by(block, prepost) %>%
+   summarise(n = n(),
+             med_samplesize = median(sample_size, na.rm = T),
+             quant25_samplesize = quantile(sample_size, .25, na.rm=T),
+             quant75_samplesize = quantile(sample_size, .75, na.rm=T),
+             IQR_samplesize = paste0(med_samplesize, " (", quant25_samplesize, ", ", quant75_samplesize, ")"),
+             med_arms = median(n_arms, na.rm = T),
+             quant25_arms = quantile(n_arms, .25),
+             quant75_arms = quantile(n_arms, .75),
+             IQR_arms = paste0(med_arms, " (", quant25_arms, ", ", quant75_arms, ")"),
+             .groups = "drop") %>%
+   select(-med_samplesize, -quant25_samplesize, -quant75_samplesize,
+          -med_arms, -quant25_arms, -quant75_arms) %>%
+   ungroup() %>%
+   mutate(pct = round((n/sum(n))*100)) %>%
+   mutate(d = paste0(pct, "\\% (", n, ")"), .after = n) %>%
+   select(-n, -pct) %>%
+   mutate(lab = case_when(
+     block == 0 & prepost == 0 ~ "Neither",
+     block == 1 & prepost == 0 ~ "Only block randomization",
+     block == 0 & prepost == 1 ~ "Only prepost",
+     block == 1 & prepost == 1 ~ "Both block randomization and prepost"
+   ), .before = block) 

> # Excluding and reorganizing columns
> design_stats <- design_stats %>%
+   select(-block, -prepost) %>%
+   select(lab, d, everything())

> tableD1 <- kbl(design_stats,
+                format = "latex",
+                align = "l",
+                escape = F,
+                col.names = c("", "Prevalence", "Median (IQR) Sample Size", "Median (IQR) Arms"),
+                booktabs = TRUE,
+                caption = "Current Use of Alternative Designs to Increase Precision\\label{tab:handcoding}") %>%
+   kable_styling(full_width = FALSE, font_size = 10) %>%
+   footnote(general = "The second column shows the percentage of designs used in a sample
+            of 217 experiments from articles published in 2022-2023, with the number of
+            experiments in parentheses. The third and fourth columns show the median sample
+            size and number of experimental arms per type of design, with the interquartile range
+            in parentheses.",
+            footnote_as_chunk = TRUE, threeparttable = TRUE)

> writeLines(tableD1, "./tables/tableD1.txt")
Runtime for 1_handcoding.R: 0.03 seconds

> #################################################################################
> # Replication file for:                                                         #
> # "Balancing Precision and Retention in Experimental Design"                    #
> #                                                                               #
> # Gustavo Diaz                                                                  #
> # Northwestern University                                                       #
> # gustavo.diaz@northwestern.edu                                                 #
> #                                                                               #
> # Erin L. Rossiter                                                              #
> # University of Notre Dame                                                      #
> # erossite@nd.edu                                                               #
> #                                                                               #
> # This file cleans the Qualtrics survey data from the three replication         #
> # experiments, produces statistics reported in Appendix E.2, and saves          #
> # the cleaned data.                                                             #
> #################################################################################
> 
> # Data ---- 
> df_dh <- readRDS("./data/raw_data/DietrichHayesReplication.rds")

> df_bg <- readRDS("./data/raw_data/BayramGrahamReplication.rds")

> df_th <- readRDS("./data/raw_data/TappinHewittReplication.rds")

> # Apply exclusion criteria -----
> 
> #  Appendix E.2 text
> 
> # Starting sample sizes
> nrow(df_dh)
[1] 1334

> nrow(df_bg)
[1] 3050

> nrow(df_th)
[1] 2598

> # Number who failed attention check
> table(df_dh$attention_check != 12 | is.na(df_dh$attention_check))

FALSE  TRUE 
 1305    29 

> table(df_bg$attention_check != 12 | is.na(df_bg$attention_check))

FALSE  TRUE 
 3024    26 

> table(df_th$attention_check != 12 | is.na(df_th$attention_check))

FALSE  TRUE 
 2570    28 

> ## Exclude those who failed attention check
> df_dh <- df_dh %>% filter(attention_check == 12)

> df_bg <- df_bg %>% filter(attention_check == 12)

> df_th <- df_th %>% filter(attention_check == 12)

> # Exclude people who did not indicate they were
> # Black or African American (value of 2)
> table(!df_dh$race %in% c(1,3,4,6,8,16)) #1.76%

FALSE  TRUE 
   23  1282 

> df_dh <- df_dh %>%
+   filter(!race %in% c(1,3,4,6,8,16))

> # Exclude people who do not have a partisan
> # affiliation or leaning from Tappin & Hewitt.
> # Article states: "(We also analyze the data excluding Independents
> # in the Appendix; the results are the same)." (pg. 54)
> table(!(df_th$pid %in% c(3,4) & df_th$pid_lean == 3)) #10.58%

FALSE  TRUE 
  272  2298 

> df_th <- df_th %>%
+   filter(!(pid %in% c(3,4) & pid_lean == 3))

> # Sample sizes
> nrow(df_dh)
[1] 1282

> nrow(df_bg)
[1] 3024

> nrow(df_th)
[1] 2298

> # Clean DH -----
> 
> df_dh <- df_dh %>%
+   # fix typo in column name
+   rename(cr_sym_whiteleg = cr_sym_whitekleg) %>%
+   # DV into one column
+   rowwise() %>%
+   mutate(speech_approve = coalesce(cr_nonsym_blackleg, cr_sym_blackleg,
+                                    cr_nonsym_whiteleg, cr_sym_whiteleg,
+                                    en_nonsym_blackleg, en_sym_blackleg,
+                                    en_nonsym_whiteleg, en_sym_whiteleg)) %>%
+   ungroup() %>%
+   # scale DV like in original study
+   # larger numbers = more approval, rescale to [0,1]
+   mutate(speech_approve_scaled = 1 - ((speech_approve - 1) / 4)) %>%
+   # clean quasi-DV same way
+   mutate(dv_quasipre_scaled = 1 - ((dv_quasipre - 1) / 4)) %>%
+   # treatment variable for comparison of interest
+   # called Z to be comparable across studies
+   mutate(Z = case_when(
+     symbolic == 1 & civilrights == 1 ~ 1,
+     symbolic == 0 & civilrights == 1 ~ 0,
+     .default = NA
+   ))

> # Clean BG -----
> 
> df_bg <- df_bg %>%
+   # DV into one column
+   rowwise() %>%
+   mutate(delpref_post = coalesce(delpref_post1, delpref_post2,
+                                  delpref_post3, delpref_post4,
+                                  delpref_post5)) %>%
+   ungroup() %>%
+   # Clean DV to 0 (US give directly), 1 (US give through IO)
+   mutate(delpref_post = case_when(
+     delpref_post == 1 ~ 0,
+     delpref_post == 2 ~ 1,
+     .default = NA
+   )) %>%
+   # Clean pre-treatment measure to match
+   mutate(delpref_pre = case_when(
+     delpref_pre == 1 ~ 0,
+     delpref_pre == 2 ~ 1,
+     .default = NA
+   )) %>%
+   # copy and rename full treatment assignment variable
+   mutate(treatment = Z) %>%
+   # treatment variable for comparison of interest
+   # called Z to be comparable across studies
+   mutate(Z = case_when(
+     treatment == "control" ~ 0,
+     treatment == "ear_cost" ~ 1,
+     .default = NA
+   ))

> # Clean TH -----
> 
> # "recode policy opinions such that higher values indicate
> # greater agreement with the in-party cue" (pg 54)
> df_th <- df_th %>%
+   # Trump disagrees. Rs need to be recoded so agreement with party cue
+   # (e.g., 1, 2...) has larger values
+   # Obama agrees. Ds do not need to be recoded. Agreement with party cue
+   # is already agreement with the question
+   mutate(across(c(salestax_pre, salestax_post, salestax_followup,
+                   pension_pre, pension_post, pension_followup,
+                   foreignaid_pre, foreignaid_post, foreignaid_followup,
+                   healthcare_pre, healthcare_post, healthcare_followup), 
+                 ~ ifelse(pid %in% 2 | pid_lean %in% 1, 8 - ., .))) %>%
+   # Obama disagrees. Ds need to be recoded so agreement with party cue
+   # (e.g., 1, 2, ...) has larger values
+   mutate(across(c(fedaudit_pre, fedaudit_post, fedaudit_followup), 
+                 ~ ifelse(pid %in% 1 | pid_lean %in% 2, 8 - ., .))) %>%
+   # make outcomes between 0 and 1 for comparability with other studies
+   mutate(pension_pre_scaled = scales::rescale(pension_pre, to = c(0, 1)),
+          pension_post_scaled = scales::rescale(pension_post, to = c(0, 1)),
+          pension_followup_scaled = scales::rescale(pension_followup, to = c(0, 1)),
+          salestax_pre_scaled = scales::rescale(salestax_pre, to = c(0, 1)),
+          salestax_post_scaled = scales::rescale(salestax_post, to = c(0, 1)),
+          salestax_followup_scaled = scales::rescale(salestax_followup, to = c(0, 1)),
+          foreignaid_pre_scaled = scales::rescale(foreignaid_pre, to = c(0, 1)),
+          foreignaid_post_scaled = scales::rescale(foreignaid_post, to = c(0, 1)),
+          foreignaid_followup_scaled = scales::rescale(foreignaid_followup, to = c(0, 1)),
+          healthcare_pre_scaled = scales::rescale(healthcare_pre, to = c(0, 1)),
+          healthcaren_post_scaled = scales::rescale(healthcare_post, to = c(0, 1)),
+          healthcare_followup_scaled = scales::rescale(healthcare_followup, to = c(0, 1)),
+          fedaudit_pre_scaled = scales::rescale(fedaudit_pre, to = c(0, 1)),
+          fedaudit_post_scaled = scales::rescale(fedaudit_post, to = c(0, 1)),
+          fedaudit_followup_scaled = scales::rescale(fedaudit_followup, to = c(0, 1))) %>%
+   # Treatment variable for comparison of interest
+   # called Z to be comparable across studies
+   mutate(Z = cue_pension, .after = cue_pension)

> # Save -----
> 
> saveRDS(df_dh, "data/processed_data/DietrichHayesReplication-clean.rds")

> saveRDS(df_bg, "data/processed_data/BayramGrahamReplication-clean.rds")

> saveRDS(df_th, "data/processed_data/TappinHewittReplication-clean.rds")
Runtime for 2_cleanreplicationdata.R: 0.42 seconds

> #################################################################################
> # Replication file for:                                                         #
> # "Balancing Precision and Retention in Experimental Design"                    #
> #                                                                               #
> # Gustavo Diaz                                                                  #
> # Northwestern University                                                       #
> # gustavo.diaz@northwestern.edu                                                 #
> #                                                                               #
> # Erin L. Rossiter                                                              #
> # University of Notre Dame                                                      #
> # erossite@nd.edu                                                               #
> #                                                                               #
> # This file conducts the block randomization for the Tappin and Hewitt          #
> # replication, saving output for use in the Wave 2 Qualtrics survey.            #
> #################################################################################
> 
> # Data -----
> 
> # Note: reading in the merged data, but when this script was originally used,
> # only the Wave 1 variables had been collected. Also, the original
> # script used the participant's Connect Identifiers. This script was adapted
> # to use the anonymous id created for the replication archive.
> df <- readRDS("./data/raw_data/TappinHewittReplication.rds")

> # Exclusions -----
> 
> ## Exclude who did not finish Wave 1
> df <- df %>% filter(finished_w1 == 1)

> ## Exclude who failed attention check
> df <- df %>% filter(attention_check == 12)

> # Clean variables for blocking -----
> 
> df <- df %>%
+   # clean demographics into dummy variables
+   mutate(gender_m = gender == 1,
+          gender_w = gender == 2,
+          gender_nb = gender == 3) %>%
+   mutate(race_white = grepl(1, race),
+          race_black_africanam = grepl(2, race),
+          race_amindian_aknative = grepl(3, race),
+          race_asian = grepl(4, race),
+          race_nativehawaiian_pi = grepl(5, race),
+          race_hispanic_latino = grepl(6, race),
+          race_arab_middleeastern = grepl(7, race),
+          race_other = grepl(8, race)) %>%
+   # include leaners in a binary partisan classification
+   mutate(pid3 = case_when(
+     pid == 1 | pid_lean == 2 ~ "D",
+     pid == 2 | pid_lean == 1 ~ "R",
+     (pid == 3 | pid == 4) & pid_lean == 3 ~ "I")
+   ) %>%
+   mutate(pid7 = case_when(
+     pid == 1 & pid_strength == 1 ~ 1, #strong Dem
+     pid == 1 & pid_strength == 2 ~ 2, #weak Dem 
+     pid_lean == 2 ~ 3, #lean Dem    
+     pid_lean == 3 ~ 4, #pure Ind
+     pid_lean == 1 ~ 5, #lean Rep  
+     pid == 2 & pid_strength == 2 ~ 6, #lean Rep    
+     pid == 2 & pid_strength == 1 ~ 7, #strong Rep    
+   )) 

> block_covars <- 
+    c(# demographics (treating age and income as continuous)
+      "age",
+      "gender_m",
+      #"gender_w", #omit bc repeats info as gender_m var
+      #"gender_nb",#omit bc no one is in design 3
+      "race_white",
+      "race_black_africanam",
+      "race_amindian_aknative",
+      "race_asian",
+      "race_nativehawaiian_pi",
+      "race_hispanic_latino",
+      "race_arab_middleeastern", 
+      "race_other",
+      "inc",
+     
+     # pre-treatment political variables
+      "pid7",
+      "ideo",
+      "ft_elites_pre_1", # Trump FT
+      "ft_elites_pre_2", # Obama FT
+      "ft_elites_pre_3", # Biden FT
+      "ft_elites_pre_4", # Romney FT
+      "ft_voters_pre_1", # Republican voters
+      "ft_voters_pre_2", # Democratic voters
+     
+     # pre-treatment measures of the outcome
+     "salestax_pre",
+     "pension_pre", # outcome of interes
+     "fedaudit_pre",
+     "foreignaid_pre",
+     "healthcare_pre")

> df_block <- df %>%
+   filter(design == 3) %>%
+   select(id, all_of(block_covars))

> # Check completeness
> # 58 incompletes excluded *prior to* block randomized treatment
> # (mentioned in Appendix Section E.3)
> completes <- complete.cases(df_block)

> table(completes)
completes
FALSE  TRUE 
   58   794 

> df_block <- df_block[completes, ]

> # Create blocks -----
> 
> # Not doing weighting or restricting the blocking wrt the main outcome
> # of interest (pension_pre) in order to use the optimal
> # algorithm (as opposed to optGreedy.)
> 
> set.seed(1234)

> out <- block(data = df_block,
+              n.tr = 2,
+              id.vars = "id",
+              block.vars = colnames(df_block)[-c(1)],
+              algorithm = "optimal",
+              distance = "mahalanobis")

> # Assign treatment -----
> 
> ## Treatment assignment within block for whether they are shown a
> ## party cue or not on each issue
> cue_salestax <- assignment(out, seed = 1, namesCol = c("cue", "nocue"))

> cue_pension <- assignment(out, seed = 2, namesCol = c("cue", "nocue"))

> cue_fedaudit <- assignment(out, seed = 3, namesCol = c("cue", "nocue"))

> cue_foreignaid <- assignment(out, seed = 4, namesCol = c("cue", "nocue"))

> cue_healthcare <- assignment(out, seed = 5, namesCol = c("cue", "nocue"))

> ## Then, treatment assignment within each block for whether
> ## they state their opinion right away or not
> state_pref <- assignment(out, seed = 6, namesCol = c("state", "wait"))

> # Reorganize -----
> 
> # Create a block id variable
> block_ids <- out$blocks$`1` %>%
+   mutate(block_id = 1:n()) %>%
+   pivot_longer(cols = `Unit 1`:`Unit 2`) %>%
+   rename(id = value)

> # Add assignments and block id to the data
> df_out <- df_block %>%
+   mutate(cue_salestax = extract_conditions(cue_salestax, df_block, "id")-1) %>%
+   mutate(cue_pension = extract_conditions(cue_pension, df_block, "id")-1) %>%
+   mutate(cue_fedaudit = extract_conditions(cue_fedaudit, df_block, "id")-1) %>%
+   mutate(cue_foreignaid = extract_conditions(cue_foreignaid, df_block, "id")-1) %>%
+   mutate(cue_healthcare = extract_conditions(cue_healthcare, df_block, "id")-1) %>%
+   mutate(state_pref = extract_conditions(state_pref, df_block, "id")-1) %>%
+   mutate(design = 3) %>%
+   plyr::join(block_ids) %>%
+   # removed unneeded blocking covariates
+   select(-(age:healthcare_pre)) %>%
+   filter(!is.na(Distance)) %>%
+   # add all ids of people to invite back from designs 1 and 2
+   bind_rows(df[,c("id", "design")] %>%
+               filter(design %in% c(1,2))) %>%
+   select(-Distance, -name)

> # Save objects -----
> 
> write.csv(df_out, file = "data/processed_data/TappinHewitt-blockrandtreatments.csv")

> # embedded data to use in Qualtrics
> # this structure: [{"WorkerID": "1234", "design ": 1},{"WorkerID": "5678", "design": 2}, ...]
> x <- paste0('{"WorkerID": "', df_out$participantId, 
+             '", "design": "', df_out$design,
+             '", "cue_salestax": "', df_out$cue_salestax,
+             '", "cue_pension": "', df_out$cue_pension,
+             '", "cue_fedaudit": "', df_out$cue_fedaudit,
+             '", "cue_foreignaid": "', df_out$cue_foreignaid,
+             '", "cue_healthcare": "', df_out$cue_healthcare,
+             '", "state_pref": "', df_out$state_pref,
+             '"',
+             '}')

> x <- paste0(x, collapse = ",")

> x <- paste0('[', x, ']')

> write_lines(x = x, file = "data/processed_data/TappinHewitt-embdata.txt")
wrote 1.00TB in  0s, 83.89PB/s                                                                                                                         Runtime for 3_tappinhewitt_blockrand.R: 0.71 seconds

> #################################################################################
> # Replication file for:                                                         #
> # "Balancing Precision and Retention in Experimental Design"                    #
> #                                                                               #
> # Gustavo Diaz                                                                  #
> # Northwestern University                                                       #
> # gustavo.diaz@northwestern.edu                                                 #
> #                                                                               #
> # Erin L. Rossiter                                                              #
> # University of Notre Dame                                                      #
> # erossite@nd.edu                                                               #
> #                                                                               #
> # This file produces all tables and figures reported in Appendix Sections       #
> # E.3, E.4, E.5, and E.6                                                        #
> #################################################################################
> 
> # Data ---- 
> df_dh <- readRDS("./data/processed_data/DietrichHayesReplication-clean.rds")

> df_bg <- readRDS("./data/processed_data/BayramGrahamReplication-clean.rds")

> df_th <- readRDS("./data/processed_data/TappinHewittReplication-clean.rds")

> # Explicit sample loss (Appendix Table E2) -----
> 
> ## Design-level sample sizes (column 3)
> ## Number of units randomized to a design, whether they finished or not
> table(df_dh$design, useNA = "ifany")

  1   2   3 
426 427 429 

> table(df_bg$design, useNA = "ifany")

   1    2    3 
1008 1010 1006 

> table(df_th$design, useNA = "ifany")

  1   2   3 
751 780 767 

> ## Lost unit at any point (column 4)
> 
> ### Dietrich & Hayes
> table(df_dh$Finished, df_dh$design, useNA = "ifany")
   
      1   2   3
  0   2   3   3
  1 424 424 426

> round(prop.table(table(df_dh$Finished, df_dh$design, useNA = "ifany"), margin = 2), 3)
   
        1     2     3
  0 0.005 0.007 0.007
  1 0.995 0.993 0.993

> ### Bayram and Graham
> table(df_bg$Finished, df_bg$design, useNA = "ifany")
   
       1    2    3
  0    5    5    7
  1 1003 1005  999

> round(prop.table(table(df_bg$Finished, df_bg$design, useNA = "ifany"), margin = 2), 3)
   
        1     2     3
  0 0.005 0.005 0.007
  1 0.995 0.995 0.993

> ### Tappin & Hewitt
> ### This study is muli-wave, so different.
> ### TRUE means they attrited somewhere *during* wave 1, 2, or 3
> ###      or they did not return for wave 2 or 3
> df_th$any_loss <- ((df_th$finished_w1 != 1) | 
+                    (df_th$finished_w2 != 1) | 
+                    (df_th$finished_w3 != 1) | 
+                    is.na(df_th$start_w2) | 
+                    is.na(df_th$start_w3))

> th_loss <- table(df_th$any_loss,
+                  df_th$design,
+                  useNA = "ifany")

> th_loss
       
          1   2   3
  FALSE 653 669 608
  TRUE   98 111 159

> round(prop.table(th_loss, margin = 2), 3)
       
            1     2     3
  FALSE 0.870 0.858 0.793
  TRUE  0.130 0.142 0.207

> ## Specifically post-treatment loss (column 5)
> 
> ### Dietrich & Hayes
> ### Because the treatment embedded data was created when the treatment
> ### screen was shown, if this variable (arbitrarily chosen 1 of 3 treatment indicators)
> ### is NA, then they attrited *prior* to treatment, if not,
> ### they attrited post-treatment
> tab_dh <- table(df_dh$design, !is.na(df_dh$civilrights), useNA = "ifany")

> tab_dh
   
    FALSE TRUE
  1     0  426
  2     2  425
  3     1  428

> round(prop.table(tab_dh, margin = 1), 3)
   
    FALSE  TRUE
  1 0.000 1.000
  2 0.005 0.995
  3 0.002 0.998

> ### Bayram and Graham
> ### Same as DH
> tab_bg <- table(df_bg$design, !is.na(df_bg$treatment), useNA = "ifany")

> round(prop.table(tab_bg, margin = 1), 3)
   
    FALSE  TRUE
  1 0.001 0.999
  2 0.002 0.998
  3 0.005 0.995

> tab_bg
   
    FALSE TRUE
  1     1 1007
  2     2 1008
  3     5 1001

> ### Tappin & Hewitt
> ### Study is muli-wave, so different.
> ### Post-treatment attrition is:
> ### (1) anytime during wave 2 for designs 1-2
> ### because treatment was immediately assigned when survey started
> ### (2) anytime during wave 2 or not starting wave 2 for design 3
> ### because treatment for design 3 was assigned between w1 and w2
> ### (3) anytime during wave 3
> ### (4) or not returning to wave 3
> 
> ### However, some participants were excluded from blocking due to 
> ### missingness in pre-treatment data.
> ### These participants definitely count as explicit sample loss,
> ### but not post-treatment loss. They were not invited back
> ### to the experimental portion or randomized to treatment.
> df_th$pt_attrit <- case_when(df_th$excluded_from_blocking == 1 ~ 0,
+                              df_th$finished_w2 == FALSE ~ 1,
+                              (df_th$design == 3) & is.na(df_th$start_w2) ~ 1,
+                              df_th$finished_w3 == FALSE ~ 1,
+                              is.na(df_th$start_w3) ~ 1,
+                              .default = 0)    

> tab_th <- table(df_th$design, df_th$pt_attrit, useNA = "ifany")

> tab_th <- tab_th[,c(2,1)]

> tab_th
   
      1   0
  1  94 657
  2 110 670
  3 101 666

> round(prop.table(tab_th, margin = 1), 3)
   
        1     0
  1 0.125 0.875
  2 0.141 0.859
  3 0.132 0.868

> # pvalues (column 6)
> fisher.test(tab_dh[c(1:2),]) #design 1 vs. 2

	Fisher's Exact Test for Count Data

data:  tab_dh[c(1:2), ]
p-value = 0.4994
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.000000 5.335101
sample estimates:
odds ratio 
         0 


> fisher.test(tab_dh[c(1,3),]) #design 1 vs. 3

	Fisher's Exact Test for Count Data

data:  tab_dh[c(1, 3), ]
p-value = 1
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
  0.00000 39.27521
sample estimates:
odds ratio 
         0 


> fisher.test(tab_bg[c(1:2),]) #design 1 vs. 2

	Fisher's Exact Test for Count Data

data:  tab_bg[c(1:2), ]
p-value = 1
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.008480817 9.631909943
sample estimates:
odds ratio 
 0.5006645 


> fisher.test(tab_bg[c(1,3),]) #design 1 vs. 3

	Fisher's Exact Test for Count Data

data:  tab_bg[c(1, 3), ]
p-value = 0.1242
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.004203937 1.782300881
sample estimates:
odds ratio 
 0.1989487 


> fisher.test(tab_th[c(1:2),]) #design 1 vs. 2

	Fisher's Exact Test for Count Data

data:  tab_th[c(1:2), ]
p-value = 0.3678
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.6407532 1.1836918
sample estimates:
odds ratio 
 0.8715334 


> fisher.test(tab_th[c(1,3),]) #design 1 vs. 3

	Fisher's Exact Test for Count Data

data:  tab_th[c(1, 3), ]
p-value = 0.759
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.6899744 1.2892153
sample estimates:
odds ratio 
 0.9434526 


> tableE2_df <- data.frame(Study = c("Dietrich \\& Hayes", "", "",
+                                  "Bayram \\& Graham", "", "",
+                                  "Tappin \\& Hewitt", "", ""),
+                        Design = c("1. Standard", "2. Pre-post", "3. Pre-post \\& blocking",
+                                   "1. Standard", "2. Pre-post", "3. Pre-post \\& blocking",
+                                   "1. Standard", "2. Pre-post", "3. Pre-post \\& blocking"),
+                        n = c(426, 427, 429,
+                              1008, 1010, 1006,
+                              751, 780, 767),
+                        losses = c("0.005 (2)", "0.007 (3)", "0.007 (3)",
+                                   "0.005 (5)", "0.005 (5)", "0.007 (7)",
+                                   "0.130 (98)", "0.142 (111)", "0.207 (159)"),
+                        pt_loss = c("0.000 (0)", "0.005 (2)", "0.002 (1)",
+                                    "0.001 (1)", "0.002 (2)", "0.005 (5)",
+                                    "0.125 (94)", "0.141 (110)", "0.132 (101)"),
+                        p = c("", "0.499", "1",
+                              "", "1", "0.124",
+                              "", "0.368", "0.759"))

> tableE2 <- kbl(tableE2_df,
+     format = "latex",
+     align = "l",
+     escape = F,
+     col.names = c("Study", "Design", "Starting sample size", 'All loss', 'Post-treatment loss', "p-value"),
+     booktabs = TRUE,
+     linesep = c("","","\\addlinespace"), #for every 3rd row space
+     caption = "Explicit Sample Loss by Study and Alternative Design  \\label{tab:exploss}") %>%
+   kable_styling(full_width = FALSE, font_size = 10) %>%
+   footnote(general = "Descriptive statistics for explicit sample loss across replication studies.
+            The starting sample size column shows the number of observations randomized to
+            each alternative design. The next two columns show the proportion and count
+            in parentheses of how many observations were lost in total and specifically
+            post-treatment. The final column shows p-values from Fisher's exact tests of
+            the null hypothesis of no difference in proportion comparing post-treatment
+            sample loss between each study's standard design and the two alternative designs.",
+            footnote_as_chunk = TRUE, threeparttable = TRUE)

> writeLines(tableE2, "./tables/tableE2.txt")

> # Prep for assessing differential sample inclusion and post-treatment attrition -----
> 
> ## clean demographics into dummy variables
> df_th_mods <- df_th %>%
+   mutate(age18_24 = age == 2,
+          age25_34 = age == 3,
+          age35_44 = age == 4,
+          age45_54 = age == 5,
+          age55_64 = age == 6,
+          age65_74 = age == 7,
+          age75_84 = age == 8,
+          age85plus = age == 9) %>%
+   mutate(gender_m = gender == 1,
+          gender_w = gender == 2,
+          gender_nb = gender == 3) %>%
+   mutate(race_white = grepl(1, race),
+          race_black_africanam = grepl(2, race),
+          race_amindian_aknative = grepl(3, race),
+          race_asian = grepl(4, race),
+          race_nativehawaiian_pi = grepl(5, race),
+          race_hispanic_latino = grepl(6, race),
+          race_arab_middleeastern = grepl(7, race),
+          race_other = grepl(8, race)) %>%
+   mutate(inc_less20 = inc == 1,
+          inc_20_34 = inc == 2,
+          inc_35_49 = inc == 3,
+          inc_50_74 = inc == 4,
+          inc_75_99 = inc == 5,
+          inc_100plus = inc == 6) %>%
+   mutate(pid_strongdem = ifelse(pid == 1 & pid_strength == 1, 1, 0),
+          pid_weakdem = ifelse(pid == 1 & pid_strength == 2, 1, 0),
+          pid_leandem = ifelse(pid_lean == 2 & is.na(pid_strength), 1, 0),
+          #pid_ind = ifelse(pid_lean == 3, 1, 0),
+          pid_leanrep = ifelse(pid_lean == 1 & is.na(pid_strength), 1, 0),
+          pid_weakrep = ifelse(pid == 2 & pid_strength == 2, 1, 0),
+          pid_strongrep = ifelse(pid == 2 & pid_strength == 1, 1, 0)) %>%
+   # political variables asked only in designs 2 and 3
+   mutate(ideo_verylib = ifelse(ideo == 1, 1, 0),
+          ideo_lib = ifelse(ideo == 2, 1, 0),
+          ideo_somewhatlib = ifelse(ideo == 3, 1, 0),
+          ideo_moderate = ifelse(ideo == 4, 1, 0),
+          ideo_somwhatcons = ifelse(ideo == 5, 1, 0),
+          ideo_cons = ifelse(ideo == 6, 1, 0),
+          ideo_verycons = ifelse(ideo == 7, 1, 0)) %>%
+   rename(ft_voters_reps = ft_voters_pre_1,
+          ft_voters_dems = ft_voters_pre_2,
+          ft_trump = ft_elites_pre_1,
+          ft_obama = ft_elites_pre_2,
+          ft_biden = ft_elites_pre_3,
+          ft_romney = ft_elites_pre_4)

> ## pretreatment demographics
> ## asked in all designs
> pretreat_demog <- df_th_mods %>% 
+   select(age18_24:pid_strongrep) %>%
+   ## remove variables with little to no variation
+   select(-age85plus, -gender_nb, -race_other,
+          -race_amindian_aknative, -race_nativehawaiian_pi,
+          -race_arab_middleeastern) %>%
+   colnames()

> ## pretreatment outcomes & extra political variables
> ## asked in designs 2 and 3
> pretreat_pol <- df_th_mods %>%
+   select(ideo_verylib:ideo_verycons,
+          ft_voters_reps:ft_romney, 
+          salestax_pre,
+          pension_pre,
+          fedaudit_pre,
+          foreignaid_pre,
+          healthcare_pre) %>%
+   colnames()

> # Differential pre-treatment sample inclusion across designs T&H (Appendix Tables E3 and E4) -----
> 
> formulas_demog <- lapply(pretreat_demog, function(var) {
+   as.formula(paste0(var, " ~ as.factor(design)"))
+ })

> formulas_pol <- lapply(pretreat_pol, function(var) {
+   as.formula(paste0(var, " ~ as.factor(design)"))
+ })

> mods_demog <- map(formulas_demog, ~ lm(.x, data = df_th_mods[!df_th_mods$any_loss,]))

> mods_pol <- map(formulas_pol, ~ lm(.x, data = df_th_mods[!df_th_mods$any_loss & df_th_mods$design %in% c(2,3),]))

> # Because there are too many models, splitting the tables
> # into two each, then combining manually as separate rows in the SI
> tableE3_top <- stargazer::stargazer(mods_demog[1:13],
+                      title = "Assessing Differential Inclusion in the Sample Across Designs - Pre-treatment Demographics Asked in All Designs",
+                      label = "tab:diffincldemog",
+                      dep.var.caption = "",
+                      dep.var.labels.include = F,
+                      star.cutoffs = .05,
+                      column.sep.width = "0pt",
+                      star.char = "*",
+                      no.space = T,
+                      font.size = "tiny",
+                      omit.stat = c("rsq", "adj.rsq", "ser", "f"),
+                      covariate.labels = c("Pre-post", "Pre-post \\& blocking", "Intercept"),
+                      column.labels =  c("Age (18-24)", "Age (25-34)", "Age (35-44)", "Age (45-54)", "Age (55-64)",
+                                        "Age (65-74)", "Age (75-84)", "Gender (Man)", "Gender (Woman)",
+                                        "Race (White)", "Race (Black or African American)",
+                                        "Race (Asian)", "Race (Hispanic or Latino)"))

% Table created by stargazer v.5.2.3 by Marek Hlavac, Social Policy Institute. E-mail: marek.hlavac at gmail.com
% Date and time: Sun, Jun 08, 2025 - 09:13:34
\begin{table}[!htbp] \centering 
  \caption{Assessing Differential Inclusion in the Sample Across Designs - Pre-treatment Demographics Asked in All Designs} 
  \label{tab:diffincldemog} 
\tiny 
\begin{tabular}{@{\extracolsep{0pt}}lccccccccccccc} 
\\[-1.8ex]\hline 
\hline \\[-1.8ex] 
 & Age (18-24) & Age (25-34) & Age (35-44) & Age (45-54) & Age (55-64) & Age (65-74) & Age (75-84) & Gender (Man) & Gender (Woman) & Race (White) & Race (Black or African American) & Race (Asian) & Race (Hispanic or Latino) \\ 
\\[-1.8ex] & (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) & (9) & (10) & (11) & (12) & (13)\\ 
\hline \\[-1.8ex] 
 Pre-post & $-$0.007 & $-$0.005 & $-$0.011 & 0.004 & 0.013 & 0.003 & 0.003 & 0.044 & $-$0.044 & $-$0.021 & 0.015 & 0.001 & 0.007 \\ 
  & (0.017) & (0.025) & (0.025) & (0.020) & (0.016) & (0.011) & (0.004) & (0.028) & (0.028) & (0.024) & (0.019) & (0.014) & (0.018) \\ 
  Pre-post \& blocking & $-$0.007 & $-$0.022 & 0.035 & $-$0.018 & 0.005 & 0.005 & 0.002 & 0.038 & $-$0.035 & $-$0.005 & 0.014 & $-$0.021 & $-$0.014 \\ 
  & (0.017) & (0.026) & (0.025) & (0.021) & (0.017) & (0.012) & (0.004) & (0.028) & (0.028) & (0.024) & (0.019) & (0.014) & (0.018) \\ 
  Intercept & 0.107$^{*}$ & 0.320$^{*}$ & 0.273$^{*}$ & 0.162$^{*}$ & 0.092$^{*}$ & 0.043$^{*}$ & 0.003 & 0.472$^{*}$ & 0.525$^{*}$ & 0.757$^{*}$ & 0.129$^{*}$ & 0.075$^{*}$ & 0.124$^{*}$ \\ 
  & (0.012) & (0.018) & (0.018) & (0.014) & (0.012) & (0.008) & (0.003) & (0.020) & (0.020) & (0.017) & (0.014) & (0.010) & (0.013) \\ 
 \hline \\[-1.8ex] 
Observations & 1,930 & 1,930 & 1,930 & 1,930 & 1,930 & 1,930 & 1,930 & 1,930 & 1,930 & 1,930 & 1,930 & 1,930 & 1,930 \\ 
\hline 
\hline \\[-1.8ex] 
\textit{Note:}  & \multicolumn{13}{r}{$^{*}$p$<$0.05; $^{**}$p$<$[0.**]; $^{***}$p$<$[0.***]} \\ 
\end{tabular} 
\end{table} 

> writeLines(tableE3_top, "./tables/tableE3_tophalf.txt")

> # (second half of table above - use only the tabular part)
> tableE3_bottom <- stargazer::stargazer(mods_demog[14:25],
+                      title = "",
+                      label = "",
+                      dep.var.caption = "",
+                      dep.var.labels.include = F,
+                      star.cutoffs = .05,
+                      column.sep.width = "0pt",
+                      star.char = "*",
+                      no.space = T,
+                      font.size = "tiny",
+                      omit.stat = c("rsq", "adj.rsq", "ser", "f"),
+                      covariate.labels = c("Pre-post", "Pre-post \\& blocking", "Intercept"),
+                      column.labels =  c("Inc (<20K)", "Inc (30-34K)", "Inc (35-49K)", "Inc (50-74K)",
+                                         "Inc (75-99K)", "Inc (100K+)", "Party (Strong Dem)",
+                                         "Party (Weak Dem)", "Party (Lean Dem)",
+                                         "Party (Lean Rep)", "Party (Weak Rep)",
+                                         "Party (Strong Rep)"))

% Table created by stargazer v.5.2.3 by Marek Hlavac, Social Policy Institute. E-mail: marek.hlavac at gmail.com
% Date and time: Sun, Jun 08, 2025 - 09:13:34
\begin{table}[!htbp] \centering 
  \caption{} 
  \label{} 
\tiny 
\begin{tabular}{@{\extracolsep{0pt}}lcccccccccccc} 
\\[-1.8ex]\hline 
\hline \\[-1.8ex] 
 & Inc (<20K) & Inc (30-34K) & Inc (35-49K) & Inc (50-74K) & Inc (75-99K) & Inc (100K+) & Party (Strong Dem) & Party (Weak Dem) & Party (Lean Dem) & Party (Lean Rep) & Party (Weak Rep) & Party (Strong Rep) \\ 
\\[-1.8ex] & (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) & (9) & (10) & (11) & (12)\\ 
\hline \\[-1.8ex] 
 Pre-post & 0.002 & 0.011 & 0.015 & $-$0.010 & 0.024 & $-$0.042 & 0.017 & 0.023 & $-$0.019 & $-$0.017 & $-$0.005 & 0.001 \\ 
  & (0.016) & (0.018) & (0.019) & (0.023) & (0.020) & (0.024) & (0.025) & (0.024) & (0.019) & (0.015) & (0.019) & (0.016) \\ 
  Pre-post \& blocking & $-$0.022 & $-$0.013 & 0.024 & 0.006 & 0.029 & $-$0.024 & $-$0.004 & 0.023 & $-$0.015 & $-$0.031$^{*}$ & 0.020 & 0.007 \\ 
  & (0.016) & (0.018) & (0.019) & (0.024) & (0.021) & (0.025) & (0.025) & (0.025) & (0.020) & (0.015) & (0.020) & (0.016) \\ 
  Intercept & 0.098$^{*}$ & 0.116$^{*}$ & 0.126$^{*}$ & 0.240$^{*}$ & 0.144$^{*}$ & 0.276$^{*}$ & 0.279$^{*}$ & 0.239$^{*}$ & 0.158$^{*}$ & 0.095$^{*}$ & 0.141$^{*}$ & 0.089$^{*}$ \\ 
  & (0.011) & (0.013) & (0.014) & (0.017) & (0.014) & (0.017) & (0.018) & (0.017) & (0.014) & (0.011) & (0.014) & (0.011) \\ 
 \hline \\[-1.8ex] 
Observations & 1,929 & 1,929 & 1,929 & 1,929 & 1,929 & 1,929 & 1,930 & 1,930 & 1,930 & 1,930 & 1,930 & 1,930 \\ 
\hline 
\hline \\[-1.8ex] 
\textit{Note:}  & \multicolumn{12}{r}{$^{*}$p$<$0.05; $^{**}$p$<$[0.**]; $^{***}$p$<$[0.***]} \\ 
\end{tabular} 
\end{table} 

> writeLines(tableE3_bottom, "./tables/tableE3_bottomhalf.txt")

> # manually add in row for p-value comparing designs 2 and 3
> ftest_pvals <- lapply(mods_demog, car::linearHypothesis, "as.factor(design)2 = as.factor(design)3")

> names(ftest_pvals) <- pretreat_demog

> ftest_pvals
$age18_24
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: age18_24 ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq     F Pr(>F)
1   1928 177.67                          
2   1927 177.67  1 1.026e-05 1e-04 0.9916

$age25_34
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: age25_34 ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   1928 413.78                           
2   1927 413.68  1  0.099775 0.4648 0.4955

$age35_44
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: age35_44 ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq     F  Pr(>F)  
1   1928 388.86                             
2   1927 388.19  1   0.67344 3.343 0.06764 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

$age45_54
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: age45_54 ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   1928 256.78                           
2   1927 256.64  1   0.14292 1.0731 0.3004

$age55_64
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: age55_64 ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   1928 170.46                           
2   1927 170.44  1   0.01837 0.2077 0.6486

$age65_74
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: age65_74 ~ as.factor(design)

  Res.Df   RSS Df  Sum of Sq      F Pr(>F)
1   1928 83.98                            
2   1927 83.98  1 0.00058875 0.0135 0.9075

$age75_84
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: age75_84 ~ as.factor(design)

  Res.Df    RSS Df  Sum of Sq      F Pr(>F)
1   1928 8.9555                            
2   1927 8.9552  1 0.00034774 0.0748 0.7845

$gender_m
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: gender_m ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   1928 481.76                           
2   1927 481.75  1  0.010814 0.0433 0.8353

$gender_w
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: gender_w ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq     F Pr(>F)
1   1928 481.82                          
2   1927 481.79  1  0.024757 0.099  0.753

$race_white
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: race_white ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   1928 364.04                           
2   1927 363.95  1  0.083786 0.4436 0.5055

$race_black_africanam
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: race_black_africanam ~ as.factor(design)

  Res.Df    RSS Df  Sum of Sq     F Pr(>F)
1   1928 229.97                           
2   1927 229.97  1 5.2414e-05 4e-04 0.9833

$race_asian
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: race_asian ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   1928 123.80                           
2   1927 123.64  1   0.15356 2.3932  0.122

$race_hispanic_latino
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: race_hispanic_latino ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   1928 207.14                           
2   1927 206.99  1   0.14508 1.3507 0.2453

$inc_less20
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: inc_less20 ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   1927 160.72                           
2   1926 160.53  1   0.19327 2.3188  0.128

$inc_20_34
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: inc_20_34 ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq    F Pr(>F)
1   1927 197.99                         
2   1926 197.81  1   0.17769 1.73 0.1886

$inc_35_49
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: inc_35_49 ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   1927 229.88                           
2   1926 229.85  1   0.02551 0.2138 0.6439

$inc_50_74
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: inc_50_74 ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   1927 350.83                           
2   1926 350.74  1  0.083241 0.4571 0.4991

$inc_75_99
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: inc_75_99 ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   1927 260.56                           
2   1926 260.56  1 0.0080618 0.0596 0.8072

$inc_100plus
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: inc_100plus ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   1927 364.55                           
2   1926 364.45  1   0.10441 0.5518 0.4577

$pid_strongdem
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: pid_strongdem ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   1928 391.95                           
2   1927 391.80  1   0.14442 0.7103 0.3995

$pid_weakdem
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: pid_weakdem ~ as.factor(design)

  Res.Df    RSS Df  Sum of Sq  F Pr(>F)
1   1928 365.37                        
2   1927 365.37  1 1.6191e-06  0 0.9977

$pid_leandem
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: pid_leandem ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   1928 241.38                           
2   1927 241.38  1 0.0052987 0.0423 0.8371

$pid_leanrep
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: pid_leanrep ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   1928 140.63                           
2   1927 140.57  1  0.058768 0.8056 0.3695

$pid_weakrep
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: pid_weakrep ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   1928 240.07                           
2   1927 239.86  1   0.20164 1.6199 0.2033

$pid_strongrep
Linear hypothesis test

Hypothesis:
as.factor(design)2 - as.factor(design)3 = 0

Model 1: restricted model
Model 2: pid_strongrep ~ as.factor(design)

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   1928 159.94                           
2   1927 159.93  1   0.01038 0.1251 0.7236


> unlist(lapply(ftest_pvals, function(x) round(x$`Pr(>F)`, 3)[2]))
            age18_24             age25_34             age35_44             age45_54             age55_64 
               0.992                0.495                0.068                0.300                0.649 
            age65_74             age75_84             gender_m             gender_w           race_white 
               0.907                0.784                0.835                0.753                0.505 
race_black_africanam           race_asian race_hispanic_latino           inc_less20            inc_20_34 
               0.983                0.122                0.245                0.128                0.189 
           inc_35_49            inc_50_74            inc_75_99          inc_100plus        pid_strongdem 
               0.644                0.499                0.807                0.458                0.399 
         pid_weakdem          pid_leandem          pid_leanrep          pid_weakrep        pid_strongrep 
               0.998                0.837                0.370                0.203                0.724 

> # Because there are too many models, splitting the tables
> # into two each, then combining manually as separate rows in the SI
> tableE4_top <- stargazer::stargazer(mods_pol[1:9],
+                      title = "Assessing Differential Inclusion in the Sample Across Designs - Additional Pre-treatment Political Questions Asked in Alternative Designs Only",
+                      label = "tab:diffinclpol",
+                      dep.var.caption = "",
+                      dep.var.labels.include = F,
+                      star.cutoffs = .05,
+                      column.sep.width = "0pt",
+                      star.char = "*",
+                      no.space = T,
+                      font.size = "tiny",
+                      omit.stat = c("rsq", "adj.rsq", "ser", "f"),
+                      covariate.labels = c("Pre-post", "Pre-post \\& blocking", "Intercept"),
+                      column.labels =  c("Ideology (Very Lib.)", "Ideology (Lib.)", "Ideology (Somewhat Lib.)",
+                                         "Ideology (Moderate)", "Ideology (Somewhat Cons.)",
+                                         "Ideology (Cons.)", "Ideology (Very Cons.)",
+                                         "Feeling Thermom. (Voters - Reps)", "Feeling Thermom. (Voters - Dems)"))

% Table created by stargazer v.5.2.3 by Marek Hlavac, Social Policy Institute. E-mail: marek.hlavac at gmail.com
% Date and time: Sun, Jun 08, 2025 - 09:13:34
\begin{table}[!htbp] \centering 
  \caption{Assessing Differential Inclusion in the Sample Across Designs - Additional Pre-treatment Political Questions Asked in Alternative Designs Only} 
  \label{tab:diffinclpol} 
\tiny 
\begin{tabular}{@{\extracolsep{0pt}}lccccccccc} 
\\[-1.8ex]\hline 
\hline \\[-1.8ex] 
 & Ideology (Very Lib.) & Ideology (Lib.) & Ideology (Somewhat Lib.) & Ideology (Moderate) & Ideology (Somewhat Cons.) & Ideology (Cons.) & Ideology (Very Cons.) & Feeling Thermom. (Voters - Reps) & Feeling Thermom. (Voters - Dems) \\ 
\\[-1.8ex] & (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) & (9)\\ 
\hline \\[-1.8ex] 
 Pre-post & $-$0.001 & 0.018 & $-$0.035 & $-$0.004 & 0.014 & $-$0.006 & 0.014 & 1.808 & 0.001 \\ 
  & (0.021) & (0.025) & (0.021) & (0.019) & (0.018) & (0.017) & (0.012) & (1.614) & (1.474) \\ 
  Pre-post \& blocking & 0.163$^{*}$ & 0.267$^{*}$ & 0.180$^{*}$ & 0.130$^{*}$ & 0.108$^{*}$ & 0.106$^{*}$ & 0.045$^{*}$ & 39.720$^{*}$ & 58.900$^{*}$ \\ 
  & (0.014) & (0.017) & (0.014) & (0.013) & (0.012) & (0.012) & (0.009) & (1.117) & (1.017) \\ 
 \hline \\[-1.8ex] 
Observations & 1,275 & 1,275 & 1,275 & 1,275 & 1,275 & 1,275 & 1,275 & 1,269 & 1,277 \\ 
\hline 
\hline \\[-1.8ex] 
\textit{Note:}  & \multicolumn{9}{r}{$^{*}$p$<$0.05; $^{**}$p$<$[0.**]; $^{***}$p$<$[0.***]} \\ 
\end{tabular} 
\end{table} 

> writeLines(tableE4_top, "./tables/tableE4_tophalf.txt")

> # (second half of table above - use only the tabular part)
> tableE4_bottom <- stargazer::stargazer(mods_pol[10:18],
+                      title = "",
+                      label = "",
+                      dep.var.caption = "",
+                      dep.var.labels.include = F,
+                      star.cutoffs = .05,
+                      column.sep.width = "0pt",
+                      star.char = "*",
+                      no.space = T,
+                      font.size = "tiny",
+                      omit.stat = c("rsq", "adj.rsq", "ser", "f"),
+                      covariate.labels = c("Pre-post", "Pre-post \\& blocking", "Intercept"),
+                      column.labels =  c("Feeling Thermom. (Elites - Trump)",
+                                         "Feeling Thermom. (Elites - Obama)",
+                                         "Feeling Thermom. (Elites - Biden)",
+                                         "Feeling Thermom. (Elites - Romney)",
+                                         "Salestax Pref.",
+                                         "Pension Pref.",
+                                         "Fed. Audit Pref.",
+                                         "Foreign Aid Pref.",
+                                         "Healthcare Pref."))

% Table created by stargazer v.5.2.3 by Marek Hlavac, Social Policy Institute. E-mail: marek.hlavac at gmail.com
% Date and time: Sun, Jun 08, 2025 - 09:13:34
\begin{table}[!htbp] \centering 
  \caption{} 
  \label{} 
\tiny 
\begin{tabular}{@{\extracolsep{0pt}}lccccccccc} 
\\[-1.8ex]\hline 
\hline \\[-1.8ex] 
 & Feeling Thermom. (Elites - Trump) & Feeling Thermom. (Elites - Obama) & Feeling Thermom. (Elites - Biden) & Feeling Thermom. (Elites - Romney) & Salestax Pref. & Pension Pref. & Fed. Audit Pref. & Foreign Aid Pref. & Healthcare Pref. \\ 
\\[-1.8ex] & (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) & (9)\\ 
\hline \\[-1.8ex] 
 Pre-post & 0.614 & $-$0.696 & $-$2.794 & $-$0.245 & 0.066 & $-$0.040 & $-$0.107 & 0.127 & 0.108 \\ 
  & (1.935) & (1.843) & (1.804) & (1.324) & (0.095) & (0.089) & (0.092) & (0.097) & (0.088) \\ 
  Pre-post \& blocking & 26.972$^{*}$ & 62.657$^{*}$ & 47.361$^{*}$ & 31.133$^{*}$ & 3.851$^{*}$ & 4.449$^{*}$ & 4.018$^{*}$ & 4.099$^{*}$ & 4.821$^{*}$ \\ 
  & (1.350) & (1.273) & (1.248) & (0.916) & (0.066) & (0.062) & (0.063) & (0.067) & (0.061) \\ 
 \hline \\[-1.8ex] 
Observations & 1,249 & 1,275 & 1,270 & 1,271 & 1,277 & 1,276 & 1,277 & 1,277 & 1,274 \\ 
\hline 
\hline \\[-1.8ex] 
\textit{Note:}  & \multicolumn{9}{r}{$^{*}$p$<$0.05; $^{**}$p$<$[0.**]; $^{***}$p$<$[0.***]} \\ 
\end{tabular} 
\end{table} 

> writeLines(tableE4_bottom, "./tables/tableE4_bottomhalf.txt")

> # Differential post-treatment attrition across designs T&H (Appendix Figures E1 and E2) -----
> 
> ## fit models for each pre-treatment covariate and design comparison
> ## although not identical indicators, removing indicator for woman because few
> ## gender non-conforming participants, so it is repetitive of "gender_m"
> formulas_demog <- lapply(pretreat_demog[-which(pretreat_demog == "gender_w")], function(var) {
+   as.formula(paste("pt_attrit ~ as.factor(design) *", var))
+ })

> formulas_pol <- lapply(pretreat_pol, function(var) {
+   as.formula(paste("pt_attrit ~ as.factor(design) *", var))
+ })

> mods_demog <- map(formulas_demog, ~ lm_robust(.x, data = df_th_mods))

> mods_pol <- map(formulas_pol, ~ lm_robust(.x, data = df_th_mods[df_th_mods$design %in% c(2,3), ]))

> ## clean and plot for demographic variables 
> 
> int_pvalues <- plyr::ldply(.data = mods_demog, .fun = function(x){
+   data.frame(p = x$p.value[5:6],
+              design = c("Pre-post vs.\nstandard design", "Pre-post & blocking vs.\nstandard design"),
+              variable = names(x$p)[4], # Variable of interest
+              color = if_else(x$coefficients[5:6] > 0, "Larger", "Smaller"),
+              row.names = NULL)
+ })

> ## No significant interactions -- gender appears significant in
> ## plot but it is not, reported in Section E.5.1
> int_pvalues$p[int_pvalues$variable == "gender_mTRUE"]
[1] 0.05794169 0.14639698

> # clean for plot
> int_pvalues <- int_pvalues %>%
+   mutate(design = factor(design, levels = c("Pre-post vs.\nstandard design", "Pre-post & blocking vs.\nstandard design")),
+          variable_plot = case_when(
+            variable == "age18_24TRUE" ~ "Age (18-24)",
+            variable == "age25_34TRUE" ~ "Age (25-34)",
+            variable == "age35_44TRUE" ~ "Age (35-44)",
+            variable == "age45_54TRUE" ~ "Age (45-54)",
+            variable == "age55_64TRUE" ~ "Age (55-64)",
+            variable == "age65_74TRUE" ~ "Age (65-74)",
+            variable == "age75_84TRUE" ~ "Age (74-84)",
+            variable == "gender_mTRUE" ~ "Gender (Man)",
+            variable == "race_whiteTRUE" ~ "Race/Ethnicity (White)",
+            variable == "race_black_africanamTRUE" ~ "Race/Ethnicity (Black or African American)",
+            variable == "race_asianTRUE" ~ "Race/Ethnicity (Asian)",
+            variable == "race_hispanic_latinoTRUE" ~ "Race/Ethnicity (Hispanic or Latino)",
+            variable == "inc_less20TRUE" ~ "Income (<20K)",
+            variable == "inc_20_34TRUE" ~ "Income (20-34K)",
+            variable == "inc_35_49TRUE" ~ "Income (35-50K)",
+            variable == "inc_50_74TRUE" ~ "Income (50-75K)",
+            variable == "inc_75_99TRUE" ~ "Income (75-99K)",
+            variable == "inc_100plusTRUE" ~ "Income (100K+)",
+            variable == "pid_strongdem" ~ "Partisanship (Strong Dem.)",
+            variable == "pid_weakdem" ~ "Partisanship (Weak Dem.)",
+            variable == "pid_leandem" ~ "Partisanship (Lean Dem.)",
+            variable == "pid_leanrep" ~ "Partisanship (Lean Rep.)",
+            variable == "pid_weakrep" ~ "Partisanship (Weak Rep.)",
+            variable == "pid_strongrep" ~ "Partisanship (Strong Rep.)",
+            .default = NA
+          ))

> ## sort based on average p-value per variable for plotting only
> avg_pvalues <- as_tibble(int_pvalues) %>%
+   group_by(variable_plot) %>%
+   summarize(avg_p = mean(p)) %>%
+   arrange(avg_p)

> int_pvalues <- int_pvalues %>%
+   mutate(variable_plot = factor(variable_plot, levels = avg_pvalues$variable_plot))

> ## plot
> plot_attrit <- ggplot(int_pvalues, aes(x = p, y = variable_plot, color = color)) +
+   facet_wrap(~design) +
+   geom_point(position = position_dodge(width = 0.2), size = 3) +
+   geom_vline(xintercept = 0.05, linetype = "dashed", color = "black") +
+   labs(x = "P-value for the difference in a characteristic's effect on attrition\nbetween the alternative and standard design",
+        y = "",
+        color = "Effect of characteristic on attrition is\nlarger/smaller in the alternative design") +
+   scale_color_manual(values = c("#2c7bb6","#abd9e9")) +
+   theme_minimal() +
+   theme(
+     axis.text.y = element_text(hjust = 1, size = 12),
+     panel.border = element_rect(color = "black", fill = NA),
+     text = element_text(size = 12),
+     axis.title = element_text(size = 12),
+     legend.title = element_text(size = 12),
+     legend.text = element_text(size = 12),
+     strip.text = element_text(size = 12),
+     legend.position = "bottom"
+   ) +
+   scale_x_continuous(breaks = c(0.05, 0.25, 0.5, 0.75, 1),
+                      labels = c("0.05", "0.25", "0.5", "0.75", "1.00"))

> ## save
> ggsave("figures/figE1.pdf", plot_attrit, width = 10, height = 6)

> ## clean and plot for political variables 
> 
> int_pvalues <- plyr::ldply(.data = mods_pol, .fun = function(x){
+   data.frame(p = x$p.value[4],
+              design = c("Pre-post & blocking vs.\nstandard design"),
+              variable = names(x$p)[3], # Variable of interest
+              color = if_else(x$coefficients[4] > 0, "Larger", "Smaller"),
+              row.names = NULL)
+ })

> # clean for plot
> int_pvalues <- int_pvalues %>%
+   mutate(variable_plot = case_when(
+            variable == "ideo_verylib" ~ "Ideology (Very Liberal)",
+            variable == "ideo_lib" ~ "Ideology (Liberal)",
+            variable == "ideo_somewhatlib" ~ "Ideology (Somewhat Liberal)",
+            variable == "ideo_moderate" ~ "Ideology (Moderate)",
+            variable == "ideo_somwhatcons" ~ "Ideology (Somewhat Conservative)",
+            variable == "ideo_cons" ~ "Ideology (Conservative)",
+            variable == "ideo_verycons" ~ "Ideology (Very Conservative)",
+            variable == "ft_voters_reps" ~ "Feeling Thermometer (Voters - Reps)",
+            variable == "ft_voters_dems" ~ "Feeling Thermometer (Voters - Dems)",
+            variable == "ft_trump" ~ "Feeling Thermometer (Elites - Trump)",
+            variable == "ft_obama" ~ "Feeling Thermometer (Elites - Obama)",
+            variable == "ft_biden" ~ "Feeling Thermometer (Elites - Biden)",
+            variable == "ft_romney" ~ "Feeling Thermometer (Elites - Romney)",
+            variable == "salestax_pre" ~ "Salestax Preference",
+            variable == "pension_pre" ~ "Pension Preference",
+            variable == "fedaudit_pre" ~ "Fed. Audit Preference",
+            variable == "foreignaid_pre" ~ "Foreign Aid Preference",
+            variable == "healthcare_pre" ~ "Healthcare Preference",
+            .default = NA
+          ))

> ## sort based on p-value size
> sort_pvalues <- as_tibble(int_pvalues) %>%
+   group_by(variable_plot) %>%
+   arrange(p)

> int_pvalues <- int_pvalues %>%
+   mutate(variable_plot = factor(variable_plot, levels = sort_pvalues$variable_plot))

> ## plot
> plot_attrit_pol <- ggplot(int_pvalues, aes(x = p, y = variable_plot, color = color)) +
+   geom_point(position = position_dodge(width = 0.2), size = 3) +
+   geom_vline(xintercept = 0.05, linetype = "dashed", color = "black") +
+   labs(x = "P-value for the difference in a characteristic's effect on attrition\nbetween the two alternative designs",
+        y = "",
+        color = "Effect of characteristic on attrition is larger/smaller\nin pre-post with blocking vs. pre-post alone") +
+   scale_color_manual(values = c("#2c7bb6","#abd9e9")) +
+   theme_minimal() +
+   theme(
+     axis.text.y = element_text(hjust = 1, size = 12),
+     panel.border = element_rect(color = "black", fill = NA),
+     text = element_text(size = 12),
+     axis.title = element_text(size = 12),
+     legend.title = element_text(size = 12),
+     legend.text = element_text(size = 12),
+     strip.text = element_text(size = 12),
+     legend.position = "bottom"
+   ) +
+   scale_x_continuous(breaks = c(0.05, 0.25, 0.5, 0.75, 1),
+                      labels = c("0.05", "0.25", "0.5", "0.75", "1.00"))

> ## save
> ggsave("figures/figE2.pdf", plot_attrit_pol, width = 10, height = 6)

> # Differential post-treatment attrition across treatment arms by design T&H (Appendix Figures E3 and E4) -----
> 
> ## fit models for each pre-treatment covariate and design comparison
> ## although not identical indicators, removing indicator for woman because few
> ## gender non-conforming participants, so it is repetitive of "gender_m"
> 
> ## also no variation in Age 74-84
> formulas_demog <- lapply(pretreat_demog[-which(pretreat_demog %in% c("gender_w", "age75_84"))], function(var) {
+   as.formula(paste("pt_attrit ~ as.factor(design) * Z *", var))
+ })

> formulas_pol <- lapply(pretreat_pol, function(var) {
+   as.formula(paste("pt_attrit ~ as.factor(design) * Z *", var))
+ })

> mods_demog <- map(formulas_demog, ~ lm_robust(.x, data = df_th_mods))

> mods_pol <- map(formulas_pol, ~ lm_robust(.x, data = df_th_mods[df_th_mods$design %in% c(2,3), ]))

> ## clean and plot for demographic variables 
> 
> int_pvalues <- plyr::ldply(.data = mods_demog, .fun = function(x){
+   data.frame(p = x$p.value[11:12],
+              design = c("Pre-post vs.\nstandard design", "Pre-post & blocking vs.\nstandard design"),
+              variable = names(x$p)[5], # Variable of interest
+              color = if_else(x$coefficients[11:12] > 0, "Larger", "Smaller"),
+              row.names = NULL)
+ })

> # clean for plot
> int_pvalues <- int_pvalues %>%
+   mutate(design = factor(design, levels = c("Pre-post vs.\nstandard design", "Pre-post & blocking vs.\nstandard design")),
+          variable_plot = case_when(
+            variable == "age18_24TRUE" ~ "Age (18-24)",
+            variable == "age25_34TRUE" ~ "Age (25-34)",
+            variable == "age35_44TRUE" ~ "Age (35-44)",
+            variable == "age45_54TRUE" ~ "Age (45-54)",
+            variable == "age55_64TRUE" ~ "Age (55-64)",
+            variable == "age65_74TRUE" ~ "Age (65-74)",
+            variable == "age75_84TRUE" ~ "Age (74-84)",
+            variable == "gender_mTRUE" ~ "Gender (Man)",
+            variable == "race_whiteTRUE" ~ "Race/Ethnicity (White)",
+            variable == "race_black_africanamTRUE" ~ "Race/Ethnicity (Black or African American)",
+            variable == "race_asianTRUE" ~ "Race/Ethnicity (Asian)",
+            variable == "race_hispanic_latinoTRUE" ~ "Race/Ethnicity (Hispanic or Latino)",
+            variable == "inc_less20TRUE" ~ "Income (<20K)",
+            variable == "inc_20_34TRUE" ~ "Income (20-34K)",
+            variable == "inc_35_49TRUE" ~ "Income (35-50K)",
+            variable == "inc_50_74TRUE" ~ "Income (50-75K)",
+            variable == "inc_75_99TRUE" ~ "Income (75-99K)",
+            variable == "inc_100plusTRUE" ~ "Income (100K+)",
+            variable == "pid_strongdem" ~ "Partisanship (Strong Dem.)",
+            variable == "pid_weakdem" ~ "Partisanship (Weak Dem.)",
+            variable == "pid_leandem" ~ "Partisanship (Lean Dem.)",
+            variable == "pid_leanrep" ~ "Partisanship (Lean Rep.)",
+            variable == "pid_weakrep" ~ "Partisanship (Weak Rep.)",
+            variable == "pid_strongrep" ~ "Partisanship (Strong Rep.)",
+            .default = NA
+          ))

> ## sort based on average p-value per variable for plotting only
> avg_pvalues <- as_tibble(int_pvalues) %>%
+   group_by(variable_plot) %>%
+   summarize(avg_p = mean(p)) %>%
+   arrange(avg_p)

> int_pvalues <- int_pvalues %>%
+   mutate(variable_plot = factor(variable_plot, levels = avg_pvalues$variable_plot))

> ## plot
> plot_attritZ <- ggplot(int_pvalues, aes(x = p, y = variable_plot, color = color)) +
+   facet_wrap(~design) +
+   geom_point(position = position_dodge(width = 0.2), size = 3) +
+   geom_vline(xintercept = 0.05, linetype = "dashed", color = "black") +
+   labs(x = "P-value for the difference in a characteristic's effect\non differential attrition across treatment arms\nbetween the alternative and standard design",
+        y = "",
+        color = "Effect of characteristic on attrition is\nlarger/smaller in the alternative design") +
+   scale_color_manual(values = c("#2c7bb6","#abd9e9")) +
+   theme_minimal() +
+   theme(
+     axis.text.y = element_text(hjust = 1, size = 12),
+     panel.border = element_rect(color = "black", fill = NA),
+     text = element_text(size = 12),
+     axis.title = element_text(size = 12),
+     legend.title = element_text(size = 12),
+     legend.text = element_text(size = 12),
+     strip.text = element_text(size = 12),
+     legend.position = "bottom"
+   ) +
+   scale_x_continuous(breaks = c(0.05, 0.25, 0.5, 0.75, 1),
+                      labels = c("0.05", "0.25", "0.5", "0.75", "1.00"))

> ## save
> ggsave("figures/figE3.pdf", plot_attritZ, width = 10, height = 6)

> ## clean and plot for political variables
> 
> int_pvalues <- plyr::ldply(.data = mods_pol, .fun = function(x){
+   data.frame(p = x$p.value[8],
+              design = c("Pre-post & blocking vs.\nstandard design"),
+              variable = names(x$p)[4], # Variable of interest
+              color = if_else(x$coefficients[8] > 0, "Larger", "Smaller"),
+              row.names = NULL)
+ })

> # clean for plot
> int_pvalues <- int_pvalues %>%
+   mutate(variable_plot = case_when(
+     variable == "ideo_verylib" ~ "Ideology (Very Liberal)",
+     variable == "ideo_lib" ~ "Ideology (Liberal)",
+     variable == "ideo_somewhatlib" ~ "Ideology (Somewhat Liberal)",
+     variable == "ideo_moderate" ~ "Ideology (Moderate)",
+     variable == "ideo_somwhatcons" ~ "Ideology (Somewhat Conservative)",
+     variable == "ideo_cons" ~ "Ideology (Conservative)",
+     variable == "ideo_verycons" ~ "Ideology (Very Conservative)",
+     variable == "ft_voters_reps" ~ "Feeling Thermometer (Voters - Reps)",
+     variable == "ft_voters_dems" ~ "Feeling Thermometer (Voters - Dems)",
+     variable == "ft_trump" ~ "Feeling Thermometer (Elites - Trump)",
+     variable == "ft_obama" ~ "Feeling Thermometer (Elites - Obama)",
+     variable == "ft_biden" ~ "Feeling Thermometer (Elites - Biden)",
+     variable == "ft_romney" ~ "Feeling Thermometer (Elites - Romney)",
+     variable == "salestax_pre" ~ "Salestax Preference",
+     variable == "pension_pre" ~ "Pension Preference",
+     variable == "fedaudit_pre" ~ "Fed. Audit Preference",
+     variable == "foreignaid_pre" ~ "Foreign Aid Preference",
+     variable == "healthcare_pre" ~ "Healthcare Preference",
+     .default = NA
+   ))

> ## sort based on p-value size
> sort_pvalues <- as_tibble(int_pvalues) %>%
+   group_by(variable_plot) %>%
+   arrange(p)

> int_pvalues <- int_pvalues %>%
+   mutate(variable_plot = factor(variable_plot, levels = sort_pvalues$variable_plot))

> ## plot
> plot_attritZ_pol <- ggplot(int_pvalues, aes(x = p, y = variable_plot, color = color)) +
+   geom_point(position = position_dodge(width = 0.2), size = 3) +
+   geom_vline(xintercept = 0.05, linetype = "dashed", color = "black") +
+   labs(x = "P-value for the difference in a characteristic's effect\non differential attrition across treatment arms\nbetween the two alternative designs",
+        y = "",
+        color = "Effect of characteristic on attrition is larger/smaller\nin pre-post with blocking vs. pre-post alone") +
+   scale_color_manual(values = c("#2c7bb6","#abd9e9")) +
+   theme_minimal() +
+   theme(
+     axis.text.y = element_text(hjust = 1, size = 12),
+     panel.border = element_rect(color = "black", fill = NA),
+     text = element_text(size = 12),
+     axis.title = element_text(size = 12),
+     legend.title = element_text(size = 12),
+     legend.text = element_text(size = 12),
+     strip.text = element_text(size = 12),
+     legend.position = "bottom"
+   ) +
+   scale_x_continuous(breaks = c(0.05, 0.25, 0.5, 0.75, 1),
+                      labels = c("0.05", "0.25", "0.5", "0.75", "1.00"))

> ## save
> ggsave("figures/figE4.pdf", plot_attritZ_pol, width = 10, height = 6)

> # Treatment effect estimates with full samples (Appendix Tables E5, E6, E7) -----
> 
> ## Design 1
> dh_design1 <- lm_robust(speech_approve_scaled ~ Z, data = df_dh, subset = design == 1)

> bg_design1 <- lm_robust(delpref_post ~ Z, data = df_bg, subset = design == 1)

> th_design1 <- lm_robust(pension_followup_scaled ~ Z, data = df_th, subset = design == 1)

> ## Design 2
> dh_design2 <- lm_robust(speech_approve_scaled ~ Z + dv_quasipre_scaled, data = df_dh, subset = design == 2)

> bg_design2 <- lm_robust(delpref_post ~ Z + delpref_pre, data = df_bg, subset = design == 2)

> th_design2 <- lm_robust(pension_followup_scaled ~ Z + pension_pre_scaled, data = df_th, subset = design == 2)

> ## Design 3
> dh_design3 <- lm_robust(speech_approve_scaled ~ Z, fixed_effects = ~ block, data = df_dh, subset = design == 3)

> bg_design3 <- lm_robust(delpref_post ~ Z, fixed_effects = ~ block, data = df_bg, subset = design == 3)

> th_design3 <- lm_robust(pension_followup_scaled ~ Z + pension_pre_scaled, fixed_effects = ~ block_id, data = df_th, subset = design == 3)

> ## Pooled treatment effects
> dh_pooled <- lm_robust(speech_approve_scaled ~ Z, data = df_dh)

> bg_pooled <- lm_robust(delpref_post ~ Z, data = df_bg)

> th_pooled <- lm_robust(pension_followup_scaled ~ Z, data = df_th)

> # Create tables
> 
> dh_mod_list <- list("Pooled" = dh_pooled, "Standard design" = dh_design1, "Pre-post" = dh_design2, "Block randomized" = dh_design3)

> bg_mod_list <- list("Pooled" = bg_pooled, "Standard design" = bg_design1, "Pre-post" = bg_design2, "Block randomized" = bg_design3)

> th_mod_list <- list("Pooled" = th_pooled, "Standard design" = th_design1, "Pre-post" = th_design2, "Block randomized" = th_design3)

> tableE5 <- modelsummary::modelsummary(dh_mod_list,
+                            coef_map = c("Z" = "Treatment Effect",
+                                         "dv_quasipre_scaled" = "Pre-treatment outcome measure",
+                                         "(Intercept)" = "Intercept"),
+                            output = "latex",
+                            stars = c('*' = .05),
+                            gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = 0)),
+                            add_rows = data.frame("Block fixed effects", "", "", "", "Yes"),
+                            title = "Treatment Effect Estimates By Design, Dietrich and Hayes Replication \\label{tab:dhests}",
+                            escape = F)

> tableE6 <- modelsummary::modelsummary(bg_mod_list,
+                            coef_map = c("Z" = "Treatment Effect",
+                                         "delpref_pre" = "Pre-treatment outcome measure",
+                                         "(Intercept)" = "Intercept"),
+                            output = "latex",
+                            stars = c('*' = .05),
+                            gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = 0)),
+                            add_rows = data.frame("Block fixed effects", "", "", "", "Yes"),
+                            title = "Treatment Effect Estimates By Design, Bayram and Graham Replication \\label{tab:bgests}",
+                            escape = F)

> tableE7 <- modelsummary::modelsummary(th_mod_list,
+                            coef_map = c("Z" = "Treatment Effect",
+                                         "pension_pre_scaled" = "Pre-treatment outcome measure",
+                                         "(Intercept)" = "Intercept"),
+                            output = "latex",
+                            stars = c('*' = .05),
+                            gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = 0)),
+                            add_rows = data.frame("Block fixed effects", "", "", "", "Yes"),
+                            title = "Treatment Effect Estimates By Design, Tappin and Hewitt Replication \\label{tab:thests}",
+                            escape = F)

> writeLines(as.character(tableE5), "./tables/tableE5.txt")

> writeLines(as.character(tableE6), "./tables/tableE6.txt")

> writeLines(as.character(tableE7), "./tables/tableE7.txt")
Runtime for 4_descriptives.R: 2.06 seconds

> #################################################################################
> # Replication file for:                                                         #
> # "Balancing Precision and Retention in Experimental Design"                    #
> #                                                                               #
> # Gustavo Diaz                                                                  #
> # Northwestern University                                                       #
> # gustavo.diaz@northwestern.edu                                                 #
> #                                                                               #
> # Erin L. Rossiter                                                              #
> # University of Notre Dame                                                      #
> # erossite@nd.edu                                                               #
> #                                                                               #
> # This file produces Figure 1.                                                  #
> #################################################################################
> 
> # Data ---- 
> df_dh <- readRDS("./data/processed_data/DietrichHayesReplication-clean.rds")

> df_bg <- readRDS("./data/processed_data/BayramGrahamReplication-clean.rds")

> df_th <- readRDS("./data/processed_data/TappinHewittReplication-clean.rds")

> # Correlations reported in Section 4.2 ----
> cor.test(df_dh$speech_approve_scaled[df_dh$design %in% c(2,3)],
+          df_dh$dv_quasipre_scaled[df_dh$design %in% c(2,3)],
+          use = "pairwise.complete.obs")

	Pearson's product-moment correlation

data:  df_dh$speech_approve_scaled[df_dh$design %in% c(2, 3)] and df_dh$dv_quasipre_scaled[df_dh$design %in% c(2, 3)]
t = 8.3558, df = 849, p-value = 2.637e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.2123900 0.3366272
sample estimates:
      cor 
0.2756594 


> cor.test(df_bg$delpref_post[df_dh$design %in% c(2,3)],
+          df_bg$delpref_pre[df_dh$design %in% c(2,3)],
+          use = "pairwise.complete.obs")

	Pearson's product-moment correlation

data:  df_bg$delpref_post[df_dh$design %in% c(2, 3)] and df_bg$delpref_pre[df_dh$design %in% c(2, 3)]
t = 35.928, df = 1346, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.6713520 0.7259379
sample estimates:
      cor 
0.6996645 


> cor.test(df_th$pension_followup_scaled[df_th$design %in% c(2,3)],
+          df_th$pension_pre_scaled[df_th$design %in% c(2,3)],
+          use = "pairwise.complete.obs")

	Pearson's product-moment correlation

data:  df_th$pension_followup_scaled[df_th$design %in% c(2, 3)] and df_th$pension_pre_scaled[df_th$design %in% c(2, 3)]
t = 33.765, df = 1278, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.6565587 0.7145574
sample estimates:
      cor 
0.6866491 


> # Simulation -----
> # Prep - create an indicator for sample loss (see 4_descriptives.R)
> df_dh$any_loss <- (df_dh$Finished == 0)

> df_bg$any_loss <- (df_bg$Finished == 0)

> df_th$any_loss <- ((df_th$finished_w1 != 1) | 
+                      (df_th$finished_w2 != 1) | 
+                      (df_th$finished_w3 != 1) | 
+                      is.na(df_th$start_w2) | 
+                      is.na(df_th$start_w3))

> # Prep - make variable for block id consistent across studies
> df_th <- df_th %>%
+   mutate(block = block_id)

> estimation_fun <- function(pct_loss, des, exp_data, dv, dv_pre, starting_n, exp){
+   
+   # estimate model under alternative designs
+   if(des == 2){
+     df_subset <- exp_data %>%
+       # starting subset used to estimate treatment effects
+       filter(design == 2 & !any_loss & !is.na(Z)) %>%
+       # choose random portion to omit to get a consistent starting n
+       # (we begin by wiping the slate clean of explicit loss)
+       sample_n(starting_n, replace = F) %>%
+       # now chose a random portion according to sim. pct_loss
+       sample_n(ceiling(starting_n - starting_n*pct_loss), replace = F)
+     
+     formula_str <- paste(dv, "~ Z + ", dv_pre)    
+     m <- estimatr::lm_robust(as.formula(formula_str),
+                               data = df_subset,
+                               subset = design == 2)
+   }
+   if(des == 3){
+     df_subset <- exp_data %>%
+       # starting subset used to estimate treatment effects
+       filter(design == 3 & !any_loss & !is.na(Z)) %>%
+       # choose random portion to delete to get to starting n
+       sample_n(starting_n, replace = F) %>%
+       # now chose a random portion according to sim. pct_loss
+       sample_n(ceiling(starting_n - starting_n*pct_loss), replace = F)
+     
+     if(exp == "th"){
+       formula_str <- paste(dv, "~ Z + ", dv_pre)  
+     }else{
+       formula_str <- paste(dv, "~ Z")  
+     }
+ 
+     m <- estimatr::lm_robust(as.formula(formula_str),
+                               fixed_effects = ~ block,
+                               data = df_subset, 
+                               subset = design == 3)
+   }
+ 
+   # store outputs
+   data.frame(exp = exp,
+              design = des,
+              n = m$n,
+              coef = m$coef["Z"],
+              se = m$std.error["Z"],
+              lower_ci = m$conf.low["Z"],
+              upper_ci = m$conf.high["Z"],
+              pct_loss = pct_loss)
+ }

> # run simulation
> inputs <- expand.grid(pct_loss = seq(0, .5, .05),
+                       des = c(2, 3))

> n_iters <- 1000

> set.seed(123)

> sim_iters_dh <- map_dfr(1:nrow(inputs), function(i) {
+   # Repeat estimation_fun n_iters times for each combination
+   map_dfr(1:n_iters, ~ estimation_fun(pct_loss = inputs$pct_loss[i],
+                                    des = inputs$des[i],
+                                    exp_data = df_dh,
+                                    dv = "speech_approve_scaled",
+                                    dv_pre = "dv_quasipre_scaled",
+                                    starting_n = 208,
+                                    exp = "dh"))
+ })

> sim_iters_bg <- map_dfr(1:nrow(inputs), function(i) {
+   # Repeat estimation_fun n_iters times for each combination
+   map_dfr(1:n_iters, ~ estimation_fun(pct_loss = inputs$pct_loss[i],
+                                   des = inputs$des[i],
+                                   exp_data = df_bg,
+                                   dv = "delpref_post",
+                                   dv_pre = "delpref_pre",
+                                   starting_n = 299,
+                                   exp = "bg"))
+ })

> sim_iters_th <- map_dfr(1:nrow(inputs), function(i) {
+   # Repeat estimation_fun n_iters times for each combination
+   map_dfr(1:n_iters, ~ estimation_fun(pct_loss = inputs$pct_loss[i],
+                                       des = inputs$des[i],
+                                       exp_data = df_th,
+                                       dv = "pension_followup_scaled",
+                                       dv_pre = "pension_pre_scaled",
+                                       starting_n = 608,
+                                       exp = "th"))
+ })

> # Standard design
> # simulate many iters at starting sample size to estimate average SE
> standard_sim_dh <- plyr::adply(1:n_iters, 1, function(i) {
+   df_subset <- df_dh %>%
+     # starting subset used to estimate treatment effects
+     filter(design == 1 & !any_loss & !is.na(Z)) %>%
+     # choose random portion to delete to get to starting n
+     sample_n(208, replace = F)
+   
+   m <- lm_robust(speech_approve_scaled ~ Z,
+                  data = df_subset,
+                  subset = design == 1)
+   
+   # store outputs
+   data.frame(n = m$n,
+              coef = m$coef["Z"],
+              se = m$std.error["Z"],
+              lower_ci = m$conf.low["Z"],
+              upper_ci = m$conf.high["Z"],
+              exp = "dh")
+ })

> standard_sim_bg <- plyr::adply(1:n_iters, 1, function(i) {
+   df_subset <- df_bg %>%
+     # starting subset used to estimate treatment effects
+     filter(design == 1 & !any_loss & !is.na(Z)) %>%
+     # choose random portion to delete to get to starting n
+     sample_n(299, replace = F)
+   
+   m <- lm_robust(delpref_post ~ Z,
+                  data = df_subset,
+                  subset = design == 1)
+   
+   # store outputs
+   data.frame(n = m$n,
+              coef = m$coef["Z"],
+              se = m$std.error["Z"],
+              lower_ci = m$conf.low["Z"],
+              upper_ci = m$conf.high["Z"],
+              exp = "bg")
+ })

> standard_sim_th <- plyr::adply(1:n_iters, 1, function(i) {
+   df_subset <- df_th %>%
+     # starting subset used to estimate treatment effects
+     filter(design == 1 & !any_loss & !is.na(Z)) %>%
+     # choose random portion to delete to get to starting n
+     sample_n(608, replace = F)
+   
+   m <- lm_robust(pension_followup_scaled ~ Z,
+                               data = df_subset,
+                               subset = design == 1)
+   
+   # store outputs
+   data.frame(n = m$n,
+              coef = m$coef["Z"],
+              se = m$std.error["Z"],
+              lower_ci = m$conf.low["Z"],
+              upper_ci = m$conf.high["Z"],
+              exp = "th")
+ })

> # Summaries
> sim_iters_all <- bind_rows(sim_iters_dh, sim_iters_bg, sim_iters_th)

> standard_sim_all <- bind_rows(standard_sim_dh,  standard_sim_bg, standard_sim_th)

> # summary for standard design
> standard_summary <- standard_sim_all %>%
+   group_by(exp) %>%
+   summarise(se_mean = mean(se), .groups = "drop")

> # summary for designs 2 and 3 to plot
> # clean up values for plot
> sim_summary <- sim_iters_all %>%
+   group_by(exp, pct_loss, design) %>%
+   summarise(se_mean = mean(se),
+             se_sd = sd(se),
+             se_max = max(se),
+             se_min = min(se),
+             se_975 = quantile(se, .975),
+             se_025 = quantile(se, .025),
+             .groups = "drop") %>%
+   # percent change relative to standard design = 
+   # [(alt design se - standard se) / standard se] × 100
+   mutate(pct_change = case_when(
+     exp == "dh" ~ ((se_mean - standard_summary$se_mean[standard_summary$exp == "dh"])/standard_summary$se_mean[standard_summary$exp == "dh"])*100,
+     exp == "bg" ~ ((se_mean - standard_summary$se_mean[standard_summary$exp == "bg"])/standard_summary$se_mean[standard_summary$exp == "bg"])*100,
+     exp == "th" ~ ((se_mean - standard_summary$se_mean[standard_summary$exp == "th"])/standard_summary$se_mean[standard_summary$exp == "th"])*100,
+   )) %>%
+   mutate(pct_change_025 = case_when(
+     exp == "dh" ~ ((se_025 - standard_summary$se_mean[standard_summary$exp == "dh"])/standard_summary$se_mean[standard_summary$exp == "dh"])*100,
+     exp == "bg" ~ ((se_025 - standard_summary$se_mean[standard_summary$exp == "bg"])/standard_summary$se_mean[standard_summary$exp == "bg"])*100,
+     exp == "th" ~ ((se_025 - standard_summary$se_mean[standard_summary$exp == "th"])/standard_summary$se_mean[standard_summary$exp == "th"])*100,
+   )) %>%
+   mutate(pct_change_975 = case_when(
+     exp == "dh" ~ ((se_975 - standard_summary$se_mean[standard_summary$exp == "dh"])/standard_summary$se_mean[standard_summary$exp == "dh"])*100,
+     exp == "bg" ~ ((se_975 - standard_summary$se_mean[standard_summary$exp == "bg"])/standard_summary$se_mean[standard_summary$exp == "bg"])*100,
+     exp == "th" ~ ((se_975 - standard_summary$se_mean[standard_summary$exp == "th"])/standard_summary$se_mean[standard_summary$exp == "th"])*100,
+   )) %>%
+   mutate(design = factor(design, levels = c(2,3),
+                          labels = c("Pre-post", "Block Randomized (and Pre-Post for Tappin & Hewitt)"))) %>%
+   mutate(exp = factor(exp, levels = c("dh", "bg", "th"),
+                       labels = c("Dietrich & Hayes",
+                                  "Bayram & Graham",
+                                  "Tappin & Hewitt")))

> implicit_loss <- ggplot(sim_summary, aes(x = pct_loss,
+                                                y = pct_change,
+                                                color = factor(design),
+                                                shape = factor(design))) +
+   facet_wrap(~exp, scales = "fixed", nrow = 1) +
+   geom_pointrange(aes(ymin = pct_change_025, ymax = pct_change_975), size = 0.3,
+                   position = position_dodge2(width = 0.03)) +
+   # Horizontal line at 0
+   geom_hline(aes(yintercept = 0),
+              linetype = "dashed",
+              color = "black") +
+   theme_minimal() +
+   theme(
+     panel.border = element_rect(color = "black", fill = NA),
+     panel.grid.minor.x = element_blank(),
+     text = element_text(size = 12),
+     axis.title = element_text(size = 12),
+     legend.title = element_text(size = 12),
+     legend.text = element_text(size = 12),
+     strip.text = element_text(size = 12),
+     legend.position = "bottom"
+   ) +
+   labs(x = "Implicit Sample Loss",
+        y = "Percentage Change in Standard Error\nRelative to Standard Design",
+        color = "Design") + 
+   scale_color_manual(values = c("grey50", "grey20")) +
+   scale_shape_manual(values = c(16,17)) + 
+   guides(color = guide_legend(title = "Design", override.aes = list(shape = c(16,17)),
+                               ncol = 2),
+          shape = "none") +
+   scale_y_continuous(limits = c(-50, 90),
+                          breaks = seq(-50, 90, by = 25),
+                          labels = paste0(seq(-50, 90, by = 25), "%")) +
+   scale_x_continuous(breaks = seq(0, .50, by = .10),
+                      labels = paste0(seq(0, 50, by = 10), "%"))

> ggsave("figures/fig1.pdf", implicit_loss, width = 10, height = 4)
Runtime for 5_replication_simulation.R: 373.38 seconds

> #################################################################################
> # Replication file for:                                                         #
> # "Balancing Precision and Retention in Experimental Design"                    #
> #                                                                               #
> # Gustavo Diaz                                                                  #
> # Northwestern University                                                       #
> # gustavo.diaz@northwestern.edu                                                 #
> #                                                                               #
> # Erin L. Rossiter                                                              #
> # University of Notre Dame                                                      #
> # erossite@nd.edu                                                               #
> #                                                                               #
> # This file conducts simulation study with data from six published experiments' #
> # replication archives.                                                         #
> #################################################################################
> 
> # Source "designer" function
> source("code/designer.R")

> # Data -----
> 
> # Publicly available data downloaded from Harvard Dataverse. See the readme
> # file for details.
> 
> # study 1
> s1 = read_dta("data/external_data/Galasso_et_al.dta")

> # study 2
> load("data/external_data/Manekin_Mitts.rdata")

> s2 = us_survey_wave2

> rm(us_survey_wave2)

> # study 3
> s3 = read.csv("data/external_data/Lyon.csv")

> # study 4
> s4 = read.csv("data/external_data/Goerger_et_al.csv")

> # study 5
> load("data/external_data/Curiel_et_al.rdata")

> s5 = imputed

> rm(imputed)

> # study 6
> s6 = read_dta("data/external_data/Simas.dta")

> # Subset and clean -----
> # Subset to experimental variables and covariates in pre-registration file
> 
> # Study 1
> df1 = s1 %>% 
+   filter(campaign != "Aggressive" & candidates_no == 2) %>% 
+   mutate(Z = ifelse(campaign == "Negative", 1, 0)) %>% 
+   select(Y = votebaldi, # outcome
+          Z, # treatment
+          Y0 = sub_riskaversion, # pre-treatment outcome proxy
+          # covariates
+          sub_trustdummy, sub_male, sub_age, sub_highschool, sub_south,
+          sub_bigtown, sub_left
+   ) %>%
+   drop_na # drop incomplete cases, if any

> # Remove single label
> attr(df1$sub_age, "label") <- NULL 

> # Study 2
> df2 = s2 %>% 
+   # Redo factors so that they don't have too many levels
+   mutate(
+     white = ifelse(race == "White / Caucasian", 1, 0)
+   ) %>% 
+   select(Y = police_action_required, # outcome 11 pt scale
+          Z = black, # treatment
+          # covariates
+          pol_views, white, partyID, female # DEVIATION FROM PREREG BECAUSE TOO MANY CATEGORIES
+   ) %>% 
+   drop_na # drop incomplete cases, if any

> # Study 3
> # IMPORTANT DEVIATION
> # We could not tell what category is represented by each number
> # so we just grouped them in arbitrarily
> df3 = s3 %>%
+   filter(treatment != "western" &
+            education < 999) %>% 
+   mutate(Z = ifelse(treatment == "african", 1, 0),
+          rel = ifelse(religion < 5, 1, 0), # There is a gap between 2 and 5
+          edu = ifelse(education < 4, 1, 0), # Edu is plausibly incremental so we pick the median value
+          mar = ifelse(marital == 2, 1, 0), # Maybe 2 is married and everything else is non-married
+          gen = gender - 1) %>% 
+   select(Y = ss_rights_index, # outcome
+          Z, # treatment
+          # covariates
+          african_welcome_homosexuals, rel, gen, edu, mar
+   ) %>% 
+   drop_na # drop incomplete cases, if any

> # Study 4
> df4 = s4 %>%
+   filter(treat <= 2) %>% 
+   mutate(Z = ifelse(treat == 2, 1, 0)) %>% 
+   select(Y = yes, # outcome
+          Z, # treatment
+          # covariates
+          actual_all_crimes5yr, officers_assaulted5yr,
+          tot_clr_murder5yr, tot_clr_index_violent5yr,
+          tot_clr_index_property5yr
+   ) %>% 
+   # coarsening for speed
+   mutate(actual_all_crimes5yr = as.character(cut(actual_all_crimes5yr, breaks = quantile(actual_all_crimes5yr), include.lowest = TRUE)),
+          officers_assaulted5yr = ifelse(officers_assaulted5yr > 0, "1", "0"),
+          tot_clr_murder5yr = ifelse(tot_clr_murder5yr > 0, "1", "0"),
+          tot_clr_index_violent5yr = as.character(cut(tot_clr_index_violent5yr, breaks = quantile(tot_clr_index_violent5yr), include.lowest = TRUE)),
+          tot_clr_index_property5yr = as.character(cut(tot_clr_index_property5yr, breaks = quantile(tot_clr_index_property5yr), include.lowest = TRUE))) %>%
+   drop_na # drop incomplete cases, if any

> # Study 5
> df5 = s5 %>%
+   select(Y = trust_score, # outcome
+          Z = assignment, # treatment
+          # covariates
+          conflict_score, voted2016, age, schooling, gender, blanco, negro, indigena
+   ) %>% 
+   drop_na # drop incomplete cases, if any

> # Study 6
> df6 = s6 %>%
+   select(Y = relative_traits, # outcome
+          Z = threat_treat, # treatment
+          # covariates
+          relative_fave, high_sj, sameparty
+   ) %>% 
+   drop_na # drop incomplete cases, if any

> # Designs -----
> designs_study1 <- expand_design(expand = T,
+                                 designer = designer,
+                                 x = list(df1),
+                                 N = nrow(df1),
+                                 drop = seq(0, .5, by = 0.1),
+                                 Y = "Y",
+                                 corr = seq(.25, .75, by = .25),
+                                 u_m = 0,
+                                 u_sd = .1,
+                                 Y_binary = TRUE,
+                                 block_many_vars = list(c("sub_trustdummy", "sub_male", 
+                                                     "sub_age", "sub_highschool", 
+                                                     "sub_south", "sub_bigtown",
+                                                     "sub_left")),
+                                 multiblock = "continuous")

> designs_study2 <- expand_design(expand = T,
+                                 designer = designer,
+                                 x = list(df2),
+                                 N = nrow(df2),
+                                 drop = seq(0, .5, by = 0.1),
+                                 Y = "Y",
+                                 corr = seq(.25, .75, by = .25),
+                                 u_m = 0,
+                                 u_sd = 1,
+                                 Y_binary = FALSE,
+                                 block_many_vars = list(c("white", "partyID", "female")),
+                                 multiblock = "discrete")

> designs_study3 <- expand_design(expand = T,
+                                 designer = designer,
+                                 x = list(df3),
+                                 N = nrow(df3),
+                                 drop = seq(0, .5, by = 0.1),
+                                 Y = "Y",
+                                 corr = seq(.25, .75, by = .25),
+                                 u_m = 0,
+                                 u_sd = .1,
+                                 Y_binary = FALSE,
+                                 block_many_vars = list( c("rel", "gen", "edu", "mar")),
+                                 multiblock = "discrete")

> designs_study4 <- expand_design(expand = T,
+                                 designer = designer,
+                                 x = list(df4),
+                                 N = nrow(df4),
+                                 drop = seq(0, .5, by = 0.1),
+                                 Y = "Y",
+                                 corr = seq(.25, .75, by = .25),
+                                 u_m = 0,
+                                 u_sd = .1,
+                                 Y_binary = TRUE,
+                                 block_many_vars = list( c("actual_all_crimes5yr", "officers_assaulted5yr",
+                                                           "tot_clr_murder5yr", "tot_clr_index_violent5yr",
+                                                           "tot_clr_index_property5yr")),
+                                 multiblock = "discrete")

> designs_study5 <- expand_design(expand = T,
+                                 designer = designer,
+                                 x = list(df5),
+                                 N = nrow(df5),
+                                 drop = seq(0, .5, by = 0.1),
+                                 Y = "Y",
+                                 corr = seq(.25, .75, by = .25),
+                                 u_m = 0,
+                                 u_sd = .5,
+                                 Y_binary = FALSE,
+                                 block_many_vars = list(c("voted2016", "age", "schooling",
+                                         "gender", "blanco", 
+                                         "negro", "indigena")),
+                                 multiblock = "continuous")

> designs_study6 <- expand_design(expand = T,
+                                 designer = designer,
+                                 x = list(df6),
+                                 N = nrow(df6),
+                                 drop = seq(0, .5, by = 0.1),
+                                 Y = "Y",
+                                 corr = seq(.25, .75, by = .25),
+                                 u_m = 0,
+                                 u_sd = 1,
+                                 Y_binary = FALSE,
+                                 block_many_vars = list(c("high_sj","sameparty")),
+                                 multiblock = "discrete")

> # set up our own simulation because DeclareDesign is actually slower
> combos <- expand.grid(drop = seq(0, .5, by = 0.1), corr = seq(.25, .75, by = .25))

> combos <- combos %>%
+   slice(rep(1:n(), each = 6))

> my_sim_function <- function(id, d){
+   ests <- draw_estimates(d)
+   ests <- ests %>%
+     mutate(sim_ID = id,
+            drop = combos$drop,
+            corr = combos$corr) %>%
+     select(-term, -statistic, -p.value, -conf.low,  -conf.high, -df, -outcome, -inquiry)
+   return(ests)
+ }

> sims <- 1000

> set.seed(1234)

> out_study2 <- plyr::adply(.data = 1:sims,.margins = 1, .fun =  my_sim_function, d = designs_study2, .id = NULL, .parallel = F)

> saveRDS(out_study2, file = "data/processed_data/diagnosis_study2.RDS")

> out_study5 <- plyr::adply(.data = 1:sims,.margins = 1, .fun =  my_sim_function, d = designs_study5, .id = NULL, .parallel = F)

> saveRDS(out_study5, file = "data/processed_data/diagnosis_study5.RDS")

> out_study6 <- plyr::adply(.data = 1:sims,.margins = 1, .fun =  my_sim_function, d = designs_study6, .id = NULL, .parallel = F)

> saveRDS(out_study6, file = "data/processed_data/diagnosis_study6.RDS")

> out_study3 <- plyr::adply(.data = 1:sims,.margins = 1, .fun =  my_sim_function, d = designs_study3, .id = NULL, .parallel = F)

> saveRDS(out_study3, file = "data/processed_data/diagnosis_study3.RDS")

> out_study1 <- plyr::adply(.data = 1:sims,.margins = 1, .fun =  my_sim_function, d = designs_study1, .id = NULL, .parallel = F)

> saveRDS(out_study1, file = "data/processed_data/diagnosis_study1.RDS")

> out_study4 <- plyr::adply(.data = 1:sims,.margins = 1, .fun =  my_sim_function, d = designs_study4, .id = NULL, .parallel = F)

> saveRDS(out_study4, file = "data/processed_data/diagnosis_study4.RDS")
Runtime for 6_simulation_study.R: 11670.09 seconds

> #################################################################################
> # Replication file for:                                                         #
> # "Balancing Precision and Retention in Experimental Design"                    #
> #                                                                               #
> # Gustavo Diaz                                                                  #
> # Northwestern University                                                       #
> # gustavo.diaz@northwestern.edu                                                 #
> #                                                                               #
> # Erin L. Rossiter                                                              #
> # University of Notre Dame                                                      #
> # erossite@nd.edu                                                               #
> #                                                                               #
> # This file produces Appendix Figure F8.                                        #                                          
> #################################################################################
> 
> # Setup -----
> 
> # ggplot global options
> theme_set(theme_bw(base_size = 20))

> # Data -----
> 
> study1 <- readRDS("data/processed_data/diagnosis_study1.RDS")

> study2 <- readRDS("data/processed_data/diagnosis_study2.RDS")

> study3 <- readRDS("data/processed_data/diagnosis_study3.RDS")

> study4 <- readRDS("data/processed_data/diagnosis_study4.RDS")

> study5 <- readRDS("data/processed_data/diagnosis_study5.RDS")

> study6 <- readRDS("data/processed_data/diagnosis_study6.RDS")

> sims <- bind_rows(study1 %>% mutate(study = "Study 1"),
+                   study2 %>% mutate(study = "Study 2"),
+                   study3 %>% mutate(study = "Study 3"),
+                   study4 %>% mutate(study = "Study 4"),
+                   study5 %>% mutate(study = "Study 5"),
+                   study6 %>% mutate(study = "Study 6")
+ )

> # Standard design -- use only when 'drop == 0' for comparison
> std = sims %>% 
+   filter(estimator == "Complete + Post-only" & drop == 0) %>%
+   select(study, corr, std_est = estimate, sim_ID)

> nrow(std)
[1] 18000

> # Left join the relevant standard design per simulation iteration
> test_df = left_join(sims, std, by = c("study", "corr",  "sim_ID"))

> # Prepare test statistics ----
> tests = test_df %>% 
+   filter(estimator != "Complete + Post-only") %>% 
+   group_by(study, estimator, corr, drop) %>% 
+   summarize(f = var.test(x = estimate, y = std_est)$statistic,
+             p_f = var.test(x = estimate, y = std_est)$p.value,
+             .groups = "drop")

> # Arrange factors for plot
> tests <- tests %>%
+   mutate(
+     estimator = factor(estimator, levels = c("Complete + Pre-post",
+                                              "Block one + Post-only",
+                                              "Block one + Pre-post",
+                                              "Block all + Post-only",
+                                              "Block all + Pre-post")),
+     corr = as.factor(corr))

> # F-test -----
> f_plot <- ggplot(tests) +
+   aes(x = drop, y = f, group = 1, color = corr, shape = p_f < 0.05) +
+   facet_grid(study ~ estimator) +
+   geom_hline(yintercept = 1, linetype = "dashed") +
+   geom_line(aes(group = corr)) +
+   geom_point(size = 2, position = position_dodge2(width = 0.03)) +
+   scale_shape_manual(values = c(4, 19)) +
+   scale_color_manual(values = c("#A0A0A0", "#6D6D6D", "#3D3D3D")) +
+   #scale_color_viridis_d() +
+   labs(x = "Sample loss",
+        y = "Variance ratio\n(alternative/standard)",
+        color = "Correlation between pre-treatment variable and Y",
+        shape = "p < 0.05") +
+   theme_minimal() +
+   theme(legend.position = "bottom",
+         panel.border = element_rect(color = "black", fill = NA),
+         panel.grid.minor.x = element_blank(),
+         panel.grid.minor.y = element_blank(),
+         text = element_text(size = 12),
+         legend.text = element_text(size = 10),
+         legend.title = element_text(size = 12),
+         legend.direction = "vertical") +
+   guides(color = guide_legend(order = 1, ncol = 5),
+          shape = guide_legend(order = 2, ncol = 2))

> # save
> ggsave(plot = f_plot, file = "figures/figF8.pdf", width = 10, height = 10)
Runtime for 7_simulation_analysis.R: 1.29 seconds

> #################################################################################
> # Replication file for:                                                         #
> # "Balancing Precision and Retention in Experimental Design"                    #
> #                                                                               #
> # Gustavo Diaz                                                                  #
> # Northwestern University                                                       #
> # gustavo.diaz@northwestern.edu                                                 #
> #                                                                               #
> # Erin L. Rossiter                                                              #
> # University of Notre Dame                                                      #
> # erossite@nd.edu                                                               #
> #                                                                               #
> # This file produces Appendix Figure G10.                                       #                                          
> #################################################################################
> 
> # Setup -----
> 
> # ggplot global options
> theme_set(theme_bw(base_size = 20))

> # Simulation -----
> 
> # Declare block randomized experiment with correlated sample loss
> # Compare with uncorrelated sample loss
> # Make blocking more consequential than in main paper
> # In this case tau is a vector of size two, one effect per block
> block_designer <- function(N, tau, p){
+   # Population
+   pop <- declare_population(
+     N = N,
+     x = rnorm(N),
+     b = ifelse(x > 0, 1, 0),
+     tau = ifelse(b == 1, tau, 0.4 - tau),
+     drop = ifelse(b == 1, 
+                   sample(c(0,1),
+                          N, 
+                          replace = T, 
+                          prob = c((1-p),p)),
+                   0),
+     # Reverse so that report = 1
+     R = ifelse(drop == 1, 0, 1)
+     )
+   
+   # Potential outcomes
+   pot_c <- declare_potential_outcomes(Y ~ tau * Z_c + x,
+                                       assignment_variables= "Z_c")
+   
+   pot_b <- declare_potential_outcomes(Y ~ tau * Z_b + x,
+                                       assignment_variables= "Z_b")
+   
+   # Estimands
+   estimand_c <- declare_inquiry(ATE_c = mean(Y_Z_c_1 - Y_Z_c_0))
+   
+   estimand_b <- declare_inquiry(ATE_b = mean(Y_Z_b_1 - Y_Z_b_0))
+   
+   # Assignments
+   assign_c <- declare_assignment(Z_c = complete_ra(N = N))
+   
+   assign_b <-declare_assignment(Z_b = block_ra(blocks = b))
+   
+   # Reveal outcomes
+   reveal_c <- declare_reveal(outcome_variables = Y,
+                              assignment_variables = Z_c)
+   
+   edit_c <- declare_step(Y_c = Y, handler = fabricate)
+   
+   reveal_b <- declare_reveal(outcome_variables = Y,
+                              assignment_variables = Z_b)
+   
+   edit_b <- declare_step(Y_b = Y, handler = fabricate)
+   
+   # Estimators
+   post_c <- declare_estimator(Y_c ~ Z_c, .method = lm_robust,
+                               inquiry = "ATE_c",
+                               label = "complete_ra")
+   
+   post_b <- declare_estimator(Y_b ~ Z_b, .method = lm_robust,
+                               inquiry = "ATE_b", 
+                               subset = R == 1,
+                               fixed_effects = ~b,
+                               label = "block_ra")
+   
+   # Return design
+   pop + 
+     pot_c + pot_b +
+     estimand_c + estimand_b +
+     assign_c + assign_b +
+     reveal_c + edit_c +
+     reveal_b + edit_b +
+     post_c + post_b
+ 
+ 
+ }

> # Expand design
> designs <- expand_design(
+   designer = block_designer,
+   N = 1000, # sample size
+   tau = c(0.2, 0.3), # Linear additive effect
+   p = seq(0, 0.8, by = 0.1), # sample loss rate
+   expand = TRUE
+ )

> # Simulate
> set.seed(20230424)

> sim_nonrandom = simulate_design(designs, sims = 1000)

> # Visualize -----
> 
> # Diagnose
> diag_nr = sim_nonrandom %>% 
+   group_by(p, estimator, tau) %>% 
+   summarize(
+     powerp = mean(p.value <= .05, na.rm = TRUE),
+     powerci = mean(0 > conf.high | 0 < conf.low, na.rm = TRUE),
+     meanp = mean(p.value, na.rm = TRUE),
+     coverage = mean(estimand <= conf.high & estimand >= conf.low, na.rm = TRUE),
+     mean_estimate = mean(estimate, na.rm = TRUE),
+     bias = mean(estimate - estimand, na.rm = TRUE),
+     mse = mean((estimate - estimand)^2, na.rm = TRUE),
+     rmse = sqrt(mean((estimate - estimand)^2, na.rm = TRUE)), 
+     mad = mean(abs(estimate - estimand), na.rm = TRUE),
+     sd_estimate = sd(estimate, na.rm = TRUE),
+     mean_se = mean(std.error, na.rm = TRUE),
+     mean_estimand = mean(estimand, na.rm = TRUE),
+     nestimand = length(unique(estimand)),
+     nsims = n()
+   ) %>% 
+   mutate(estimator = recode(estimator,
+                             "block_ra" = "Block",
+                             "complete_ra" = "Complete"),
+          tau = recode(tau,
+                       `0.3` = "(0.3, 0.1)",
+                       `0.2` = "(0.2, 0.2)"))

> # recode order of estimators
> diag_nr$estimator = fct_relevel(diag_nr$estimator,
+                                 "Block",
+                                 "Complete")

> diag_nr$power_label = "Power"

> diag_nr$bias_label = "Bias"

> # Plot
> power_nr = ggplot(diag_nr) +
+   aes(x = p, y = powerp, shape = as.factor(tau), color = estimator) +
+   geom_point(size = 3) + geom_line(size = 1) +
+   theme(legend.position = "bottom") +
+   labs(x = "Sample loss rate",
+        y = expression(paste('Power at ', alpha  == 0.05)),
+        color = "Randomization",
+        shape = "True effects") +
+   scale_color_grey(start = 0, end = 0.5) +
+   facet_wrap(~ power_label) +
+   scale_y_continuous(breaks = seq(0, 1, by = 0.2), limits = c(0, 1))

> bias_nr = ggplot(diag_nr) +
+   aes(x = p, y = abs(bias), shape = as.factor(tau), color = estimator) +
+   geom_point(size = 3) + geom_line(size = 1) +
+   theme(legend.position = "bottom") +
+   labs(x = "Sample loss rate",
+        y = "Aboslute bias: |Mean(estimate - estimand)|",
+        color = "Randomization",
+        shape = "True effects") +
+   scale_color_grey(start = 0, end = 0.5) +
+   facet_wrap(~ bias_label)

> power_bias = ggarrange(power_nr, bias_nr, ncol = 2, nrow = 1, 
+                        common.legend = TRUE, legend = "bottom")

> power_bias

> #save
> ggsave(plot = power_bias, filename = "figures/figG10.pdf",
+        height = 6, width = 12)
Runtime for 8_simulation_nonrandom.R: 331.83 seconds
Replication completed at:  2025-06-08 12:39:53.03735 
